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1.  Introduction 


Chemical  vapor  deposition  (CVD)  is  arguably  one  of  the  most  important  techniques  for 
processing  of  structural  materials.  The  classification  ‘structural  materials’  here  includes  high 
temperature  aerospace  materials  such  as  composite  structures,  ceramic-matrix  composites 
(CMC),  functionally-gradient  materials  (FGM),  etc.  Representative  examples  are  SiC  (fiber, 
coating),  Si3N4,  BN,  other  carbides  and  nitrides,  composites  such  as  SiC-SiC,  carbon-SiC, 
carbon-carbon,  and  carbon-BN. 

Several  variations  of  CVD  are  used  in  the  preparation  of  the  above  composite  materials. 

(1)  Conventional  CVD  on  flat  substrates:  example,  SiC,  Si3N4.  Here,  the  processing  is 
similar  to  that  in  electronics  and  optical  applications. 

(2)  CVD  coating  of  complex  arbitrary  shapes  and  internal  cavities  such  as  in  the  following 
examples.  Hard  CVD  coatings  are  used  in  rocket  nozzles  and  gas  turbine  engines  for 
erosion  resistance.  TiB2  and  TiC  coatings  are  used  for  coating  titanium  blades  and  Si3N4 
coatings  on  carbon  nozzles.  AI2O3,  TiN,  TiC  and  TiB2  coatings  are  used  on  tool  bits  to 
increase  the  cutting  life  of  tools  through  improvements  in  hardness  and  wear  resistance. 

(3)  CVD  of  fibers.  Examples  include  SiC,  Si3N4,  boron  and  oxide  fibers.  These 
monofilament  fibers  are  mainly  used  for  reinforcing  with  other  materials  to  produce 
desirable  composites. 

(4)  CVD  of  preforms.  Here  the  CVD  process  is  used  to  infiltrate  a  preform  made  from 
one  of  the  above  fibers  to  form  a  composite.  Composites  such  as  SiC-SiC,  C-C,  C-BN, 
and  C-SiC  have  been  produced  by  this  technique.  CVD  inside  preforms  is  generally  called 
chemical  vapor  infiltration  (CVI). 

In  this  work,  we  consider  all  of  the  above  four  variations  of  CVD,  including  CVI  and 
related  applications.  The  superior  thermal  stability  of  ceramic  matrix  composites  (CMC) 
compared  to  polymers  and  metals  makes  them  unique  for  many  high  temperature  applications. 
This  fact,  combined  with  low  density  and  chemical  inertness,  make  CMCs  attractive  for  many 
aerospace  applications.  Carbon-carbon  composites  are  used  in  rocket  nozzles,  nosecones  for  re¬ 
entry  vehicles,  heat  shields,  and  aircraft  brakes.  SiC-SiC  and  C-SiC  find  many  high  temperature 
applications.  Zr02  composite  reinforced  with  TiN  whiskers  is  a  good  candidate  for  advanced  heat 
engines  due  to  a  good  match  of  thermal  properties.  In  producing  the  above  composites, 
improvements  in  reliability  and  performance  can  be  expected  when  such  materials  are  designed 
with  a  continuous  compositional,  and  therefore,  fijnctional  gradient.  Further  details  on  the  above 
discussion  can  be  found  in  a  recent  monograph  by  Galasso  [1],  Savage  [2]  and  a  National 
Research  Council  report  [3]. 

Despite  the  impressive  range  of  applications  outlined  above,  a  fundamental  understanding 
of  CVD  of  structural  materials  is  lacking.  A  recent  National  Research  Council  report  [3]  on  high 
performance  composites,  prepared  by  a  distinguished  panel  for  the  National  Materials  Advisory 
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Board,  states;  “A  basic  science  approach  to  the  CVD...,  providing  fundamental  insights  into 
decomposition  kinetics,  fiber-formation  mechanisms,  and  surface  reactions,  can  result  in  a  better 
optimization  of  processes  and  products.  Although  it  may  not  be  inordinately  difficult  to  form  a 
new  CVD  fiber  of  a  given  composition  initially,  achievement  of  consistency  in  properties  may 
require  understanding  and  carefully  controlling  the  parameters  of  the  deposition  process.” 

The  above  conclusion  is  even  more  true  when  the  composition  is  varied  to  produce 
functionally  gradient  materials  [FGM].  In  addition,  process  scale-up  is  difficult  because  of  the 
simultaneous  heat,  momentum  and  mass  transfers,  along  with  gas  phase  and  surface  reactions,  and 
complex  flow-induced  orientation  effects.  In  this  regard,  mathematical  modeling  of  the  transport 
phenomena  associated  with  the  process,  and  computer  simulation  will  be  valuable  in  many 
respects,  as  listed  below: 

(1)  to  provide  an  understanding  of  the  process  mechanisms; 

(2)  to  develop  a  cause  and  effect  picture  between  the  process  variables  (for  which  there  are 
knobs  on  the  reactor  panel)  vs.  process  figures-of-merit,  such  as  growth  rate,  uniformity, 
surface  material  composition,  etc.  (note  that  these  figures-of-merit,  in  turn,  affect  the 
structural  characteristics,  including  mechanical  and  thermal  properties); 

(3)  process  optimization  to  meet  production  goals; 

(4)  equipment  design  to  meet  given  processing  objective;  and 

(5)  process  control 

That  is,  the  models  and  computer  code  can  serve  as  a  cost-effective  design  tool.  There  has  been 
very  little  modeling  work  for  CVD  of  composites.  While  it  is  true  that  CVD  modeling  has  been 
around  for  about  a  decade  (mainly  devoted  to  electronics  materials)  and  this  knowledge  base 
would  be  valuable  for  structural  materials  as  well,  straightforward  application  of  commercial 
computer  codes  to  CVD  of  aerospace  materials  is  not  possible  at  present.  There  are  several 
unique  issues  which  need  innovations  in  terms  of  both  physical  models  and  suitable  computational 
techniques.  These  include; 

•  models  for  multicomponent  mass  transfer  with  chemical  reactions  in  composites  processing, 

•  surface  processes, 

•  sound  models  for  thermal  diffusion, 

•  preform  infiltration, 

•  model  to  represent  the  random  pore  network  of  the  preforms, 

•  development  of  suitable  boundary  conditions, 

•  transient  analysis  to  study  growth  of  FGMs 

•  computational  techniques  appropriate  for  the  above  physical  processes  both  in  terms  of 
simulation  speed  and  accuracy,  and 

•  models  to  correlate  the  transport  phenomena  effects  to  the  microstructural  mechanical  and 
thermal  properties. 
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It  is  the  intention  of  the  present  STTR  work  to  develop  suitable  transport  and  reaction 
models  for  CVD  and  CVI  of  composites  and  thus  a  computational  simulation  tool.  This  work  has 
signihcant  relevance  to  Air  Force  missions  since  high  temperature  composites  are  extensively  used 
in  militaiy  applications.  Improved  processes  through  better  understanding  provided  by  modeling 
can  impact  on  the  material  quality  and  cost  of  processing.  Commercial  applications  for 
composites  are  also  impressive  in  commercial  aviation,  high  speed  civil  transport,  automotive, 
tool  industry,  and  power  generation,  to  name  a  few.  Since  the  code  and  model  development  are 
done  in  a  modular  fashion,  application  to  electronics  materials  processing  is  readily  possible  (and 
is  demonstrated  here  through  selected  simulations  of  silicon  epitaxy).  Hence,  there  is  a  market  for 

a  user-friendly  computer  code  for  CVD/CVI  simulation  and  we  intend  to  market  the  code  after 
this  STTR. 

The  results  from  the  Phase  I  STTR  project  are  presented  in  this  report.  This  final  report  is 
organized  as  follows:  Section  2  lists  the  Phase  I  technical  objectives  as  given  in  our  proposal  and 
a  summary  of  accomplishments.  Section  3  provides  a  background  to  CVD  and  CVI  processes  to 
prepare  for  the  framework  of  analysis  presented  in  Sections  4  and  5.  The  analysis  in  Section  4 
discusses  the  CVD  modeling  and  code  development;  only  those  features  which  were  incorporated 
in  our  model  and  code  in  Phase  I  are  discussed  at  great  length;  other  standard  (by  now)  features 
are  briefly  discussed  for  the  sake  of  completion.  This  is  followed  by  a  similar  exposition  for  CVT 
in  Section  5.  Section  6  is  devoted  to  the  results  of  Phase  I  work  and  discusses  applications  of 
boron  fiber  coating,  SiC  growth,  and  silicon  epitaxy;  finally,  results  for  the  CVT  of  composites  are 
presented  in  Section  7.  Concluding  remarks  and  an  overview  of  our  Phase  II  directions  are  given 
in  Section  8. 
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2.  Phase  I  Objectives  and  Accomplishments 


The  objective  of  this  multiphase,  multiyear  STTR  program  is  to  conduct  comprehensive 
transport  phenomena  and  kinetic  studies  through  computer  modeling,  in  chemical  vapor 
deposition  and  infiltration  of  composite  materials.  Mathematical  models  are  to  be  developed  for 
heat  and  momentum  transfer,  multicomponent  species  transport  with  chemical  reactions,  surface 
phenomena,  infiltration  inside  preforms,  and  correlation  between  macroscopic  phenomena  and 
microscopic  structural  properties.  These  models  are  to  be  incorporated  in  an  existing  SRA 
computational  fluid  dynamics  (CFD)  code.  A  CVD  specific  code  and  graphical  user  interface 
(GUI)  are  to  be  developed.  The  code  and  model  are  to  be  used  to  study  CVD  and  CVI 
processing  of  representative  composites.  The  theoretical  work  needs  to  be  validated  by 
comparison  with  experiments.  Specific  growth  work  in  progress  at  Air  Force  Laboratories  will  be 
studied  using  the  code  with  a  view  to  optimize  processes  and  develop  processes  for  new  material 
systems.  The  motivation  for  this  STTR  project  comes  from  the  current  need  to  develop 
computational  simulation  tools  which  can  serve  as  cost-effective  design  aids  in  development, 
optimization,  and  control  of  CVD  and  CVI  processes  in  the  preparation  of  aerospace  structural 
materials. 

The  Phase  I  scope  is  necessarily  limited  and  the  specific  Phase  I  objectives  are: 

(1)  to  develop  and  incorporate  suitable  transport  models  for  CVD  in  an  existing  SRA 
computational  code; 

(2)  to  demonstrate  the  validity  of  the  model  by  performing  simulations  of  CVD  processing  of 
SiC-based  composites;  and 

(3)  to  develop  preliminary  models  for  the  infiltration  of  porous  preforms  in  processing  ceramic 
matrix  composites  and  functionally  gradient  materials. 

The  above  Phase  I  objectives  have  been  achieved  successfiilly. 

(1)  A  sound  model  for  thermal  diffusion  coefficients  based  on  first-principles  has  been 
incorporated  and  verified. 

(2)  Boundary  conditions  for  species  which  incorporate  multicomponent  diffusion  and  thermal 
diffusion  have  been  developed  and  coded.  These  conditions  also  include  a  proper 
description  of  surface  deposition  mechanisms. 

(3)  A  formal  procedure  to  handle  surface  species  has  been  developed. 

(4)  In  structural  materials  processing,  the  feed  gases  decompose  into  smaller  fi'agments 
consisting  of  radicals  and  atoms.  It  is  not  uncommon  to  encounter  20  or  more  species 
participating  in  as  many  or  more  reactions.  Such  a  situation  may  be  overwhelming  in 
multidimensional  simulations  since  it  drives  up  the  number  of  equations  to  be  solved. 

(Note:  each  species  needs  a  mass  conservation  equation.)  We  have  adapted  an  existing 
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SRA  well-mixed  reactor  code  (zero-dimension)  and  another  SRA  one-dimensional  code 
for  the  purpose  of  developing  an  “abbreviated  chemistry  set”  for  transfer  to 
multidimensional  simulations.  That  is,  a  given  application  is  first  run  using  these  “lower- 
order”  codes  with  the  same  process  input  but  with  a  very  large  number  of  species  and 
elementary  chemical  reactions.  The  results  are  carefully  analyzed  for  identification  of 
trace  species,  and  reactions  which  contribute  only  marginally  in  the  time-scale  of  interest. 
We  ensure  that  the  overall  characteristics  do  not  change  appreciably  when  such  species 
and  reactions  are  dropped  from  the  set.  Then  a  smaller,  manageable  set  is  transferred  to 
the  multidimensional  model. 

(5)  On  the  heat  transfer  side,  we  incorporated  a  “conjugate  heat  transfer”  module  into  the 
CVD  analysis  code.  This  module  allows  examination  of  heat  transfer  within  the  solid 
(such  as  a  wall,  wire,  object  being  coated)  which  is  in  contact  with  the  fluid  stream.  Such 
a  coupled  analysis  is  self-consistent  and  needed  in  situations  where  the  solid  in  question  is 
heated  with  a  specified  heat  flux  and  the  resistance  to  transport  Avithin  that  solid  has  an 
impact  on  the  heat  transfer  in  the  fluid. 

(6)  Finally,  models  for  the  preform  infiltration  have  been  developed.  This  aspect  was  the 
responsibility  of  our  STTR  university  partner — ^the  University  of  Houston. 

(7)  The  models  and  computer  code,  including  all  of  the  above  developments,  have  been 
demonstrated  by  performing  the  following  computations:  (a)  CVD  coating  of  fibers  with 
boron,  (b)  SiC  growth  by  CVD,  and  (c)  infiltration  of  preforms  with  SiC.  Selected 
computations  of  silicon  epitaxy  were  also  performed  to  demonstrate  the  applicability  of 
our  work  in  the  electronics  materials  field;  it  is  our  understanding  that  the  latter  is  also  of 
interest  to  AFOSR  and  Air  Force  Laboratories. 

(8)  It  is  important  to  address  practical  issues  of  concern  to  the  structural  materials  community 
for  the  proposed  modeling  work  to  have  a  significant  impact.  SRA  has  initiated 
interaction  with  three  large  companies  involved  in  the  preparation  of  structural 
composites:  Pratt  &  Whitney  for  CVI  processing  of  composites;  Textron  Specialty 
Materials  on  high  strength  fiber  CVD  of  boron,  SiC  and  other  materials  and  rapid  CVI 
densification  of  several  composites;  and  Morton  Advanced  Materials,  specializing  in  SiC 
CVD  for  nozzles  and  many  aerospace  applications.  All  three  companies  have  helped  to 
define  problems  and  will  provide  data  for  model  validation  in  Phase  II. 
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3.  Background  on  CVD  and  CVI 


In  this  section,  we  provide  a  background  on  CVD  and  CVI  processes  and  equipment, 
modeling  issues,  and  modeling  status  prior  to  this  Phase  I  work. 

3.1  CVD 


CVD  is  a  process  wherein  feed  gases  react  in  the  gas  phase  to  form  precursors  which  react 
on  and  with  the  substrate  surface  to  deposit  a  solid  film  of  desirable  composition  on  the  substrate. 
The  substrate  is  usually  heated  to  high  temperatures  (>500°C).  A  variety  of  commercial  reactors 
is  available  which  vary  in  orientation  (horizontal  vs.  vertical),  single  vs.  multiwafer,  cold  wall  vs 
hot  wall,  substrate  rotation,  type  of  heating,  etc.  Most  CVD  equipment  for  electronics  materials 
and  some  structural  materials  are  designed  for  deposition  on  flat  substrates/wafers.  In  large  scale 
CVD  coating  of  complex  objects,  the  reactor  geometry  and  size,  gas  handling  system,  object 
holding  pedestal  and  its  rotation  are  significantly  different  from  those  used  in  the  semiconductor 
industry.  The  reactors  are  often  one  of  a  kind  and  built  for  the  needs  at  hand. 

In  CVD  on  fibers,  the  fiber  substrate  continuously  moves  through  the  reactor.  The  reactor 
is  usually  a  long,  vertical  or  horizontal  cylinder  made  of  quartz  or  similar  materials.  The  reactor 
diameter  is  normally  small  (<2”)  since  the  fiber  to  be  coated  is  less  than  a  millimeter  in  diameter, 
therefore  making  it  economical  in  terms  of  gas  handling,  heating,  etc.  Heating  the  fiber  may  be 
achieved  by  an  electric  current  through  the  fiber  or  induction  heating  of  the  reactor.  Typical  fiber 
temperature  for  deposition  of  high  strength  materials  can  be  in  the  range  of  1000-2000°C. 

Most  of  the  early  CVD  modeling  has  been  for  electronics  materials;  modeling  work  on 
CVD  of  composites  is  very  limited.  Much  of  the  progress  in  CVD  modeling  has  been  pioneered 
by  Jensen  in  a  series  of  studies  [4].  These  studies  involved  solution  of  Navier-Stokes  equations 
and  energy  equation  for  the  carrier  gas.  For  species  transport,  usually  mass  transfer  of  a  single 
key  species  was  assumed  to  be  important  and  a  conservation  equation  was  solved  for  this  species. 
This  approach  was  used  to  study  horizontal  and  vertical  reactors  by  Jensen  and  later  by  many 
others.  A  thorough  review  of  the  CVD  models  in  the  literature  is  provided  in  a  new  book  by  the 
Principal  Investigator  (PI)  of  the  current  study  [5].  These  models  have  provided  the  necessary 
understanding  on  the  role  of  the  boundary  layer,  mass  transfer  vs.  kinetic  control,  effect  of 
buoyancy  and  recirculating  flows,  thus  shaping  the  development  of  equipment  and  processes. 

Further  progress  is  currently  needed  in  several  areas.  The  first  issue  is  to  incorporate 
detailed  gas  phase  chemistry  in  multidimensional  transport  simulations.  Data  on  reaction 
pathways  and  rate  coefficients  are  increasingly  becoming  available  both  from  experiments  and 
quantum  chemistry  models  (though  such  development  for  structural  materials  lags  behind  that  for 
electronics  materials).  While  there  are  several  studies  incorporating  a  large  number  of  species 
reactions,  the  corresponding  transport  modeling  is  typically  one-dimensional  impinging  jet  or 
rotating  disk  geometries,  [for  example.  Ref  6  by  the  Sandia  researchers  on  silicon  CVD].  We 
need  detailed  chemistry  modeling  along  with  2-  and  3-dimensional  transport  simulations  for 
structural  materials.  Development  of  surface  chemistry  modeling  lags  even  more.  Here,  rather 
than  assuming  sticking  coefficients  for  radical  recombination  on  the  surface,  the  mechanism  must 
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be  considered  in  terms  of  elementaiy  steps  involving  species  diffusion,  adsorption,  surface  site 
creation  and  coverage,  surface  reaction,  and  desorption.  Appropriate  adsorption  isotherms  (such 
as  Langmuir-Hinshelwood  kinetics,  Eley-Ridel  kinetics)  must  be  incorporated.  A  sound  model 
for  computing  thermal  diffusion  coefficients  based  on  first  principles  is  needed  to  overcome  the 
dependence  on  empirical  relations  which  are  hard  to  obtain.  Development  of  meaningful 
boundary  conditions,  especially  for  species  transport  and  heat  transfer,  needs  attention.  For 
composite  materials,  the  effect  of  process  conditions  and  the  macroscopic  transport  phenomena 
on  microscopic  properties  such  as  thermal  stress  in  the  composite  must  be  ascertained.  Finally,  a 
thorough  model  validation  by  comparison  with  experimental  results  is  critical  to  have  confidence 
in  the  predictive  capabilities  of  the  model  for  design  purposes.  All  these  issues  are  the  focus  of 
this  work. 

In  many  of  the  above  cases,  innovations  are  needed  not  only  in  models  but  also  in 
computational  techniques;  it  is  the  lack  of  the  latter,  to  some  extent,  that  has  prevented  a 
widespread  use  of  computational  codes  for  structural  materials  applications.  In  this  regard,  we 
note  that  there  is  extensive  development  in  mathematical  analyses  and  CFD  techniques  in  related 
fields  such  as  combustion  and  aerodynamics.  These  developments  must  be  taken  advantage  of  in 
CVD  modeling.  In  this  regard,  existing  commercial  codes  such  as  FLUENT,  PHOENICS, 
FIDAP,  and  our  own  MINT  have  made  significant  strides.  Future  work  must  focus  on 
innovations  in  algorithms  to  handle  the  computational  burden  imposed  by  a  large  number  of 
species,  multicomponent  effects,  thermal  diffusion,  variable  species  properties,  transient 
simulation  capabilities  for  composite  materials,  techniques  for  efficient  3-D  simulations,  and 
graphical  user-interface  specifically  geared  toward  CVD  and  the  materials  community.  Such 
issues  and  development  work,  provided  that  the  solutions  are  innovative,  are  ideally  suited  for  an 
STTR  project  and  can  result  in  a  fruitful  technology  transfer  to  the  private  sector. 

3.2  CVI 


Next  we  discuss  issues  related  to  CVI  of  composites.  Typically  in  CVI,  reactant  gases 
diffuse  through  the  pores  of  a  preform  and  react  on  the  surfaces  of  the  fibers  which  lead  to  the 
filling  of  the  preform  with  the  matrix  material.  The  isothermal  CVI  reactor  is  somewhat  similar  to 
a  conventional  CVD  reactor  (vertical  or  horizontal)  and  contains  a  preform  supported  on  a 
susceptor.  The  preform  can  be  any  arbitrary  3-D  structure.  The  preform  structure  is  heated  by 
external  means  (induction,  rf  or  any  other  heating).  There  is  very  little  pressure  drop  in  the 
reactor  to  force  the  flow  through  the  fiber  bundle.  In  isothermal  CVI  reactors,  the  'densification' 
of  the  preform  is  usually  slow  (order  of  days)  since  diffusion  is  the  dominant  mechanism.  In 
addition,  premature  closing  of  the  pores  can  occur  resulting  in  density  gradients  (reason  explained 
below)  and  excessive  porosity.  These  features,  in  any  case,  are  function  of  temperature,  pressure, 
thermodynamics  and  kinetics  of  the  system  [1]. 

Commonly  in  many  systems,  the  reaction  rate  (gas  phase  or  surface  reaction)  increases 
with  temperature.  In  isothermal  CVI,  conventional  heating  leads  to  higher  temperatures  close  to 
the  periphery  of  the  preform  with  a  cooler  interior.  Such  a  temperature  gradient  in  the  ‘wrong 
direction’  results  in  rapid  reactions  first  on  the  outer  side  resulting  in  premature  pore  closure. 
Hence,  it  is  not  unusual  to  take  the  preform  out  of  the  reactor  after  some  time,  grind  the  surface. 
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place  it  back  in  the  reactor,  and  reimpregnate  to  obtain  good  matrix  penetration  [1].  Actually 
what  is  ideal  is  a  temperature  profile  that  will  give  ‘inside-out’  densification  to  reduce  processing 
times.  By  controlling  the  heating  mechanism  and  hence,  the  temperature  profile,  this  can  be 
achieved.  Professor  Economou,  our  STTR  university  partner,  has  shown  using  his  models  that 
microwave  heating  or  pulsed  heating  with  external  cooling  can  provide  desirable  temperature 
profiles  in  the  ‘right  direction’.  With  such  modifications,  conventional  no-flow  CVI  still  is  a 
method  of  potential  for  a  variety  of  applications  despite  long  processing  times.  Other  techniques 
to  accelerate  the  process  time  include  application  of  pressure  and  temperature  gradients  [7]. 
Applying  a  pressure  gradient  forces  a  flow  through  the  preform.  Such  gradient  processes  have 
indeed  reduced  the  processing  time,  though  at  present  these  methods  can  only  process  small 
objects. 


The  modeling  of  infiltration  inside  the  preform  has  been  considered  by  Professor 
Economou  using  the  ‘dusty  gas  model’  approach  (described  in  section  5)  and  others  using 
different  representations  for  the  porous  network  [8-11].  Critical  issues  in  preform  modeling  relate 
to  representation  of  the  random  porous  network  of  the  preform,  models  for  the  flow  through  the 
fiber  geometry  and  reactions  inside  the  pores,  and  finally,  coupling  the  preform  model  to  the 
reactor  conditions  outside  the  preform.  Further  discussion  of  CVI  modeling  issues  and  details  is 
postponed  to  Section  5. 
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4.  CVD  Reactor  Analysis 


The  analysis  described  in  this  section  applies  to  CVD  reactors  used  in  conventional 
deposition,  coating  of  objects  of  arbitrary  size  and  shape,  and  fibers.  It  also  applies  to  CVI 
reactors  for  modeling  conditions  external  to  the  preform;  some  exceptions  and  specific  models  for 
infiltration  into  the  preform  are  subjects  of  Section  5.  For  details  beyond  the  discussion  given 
below,  reference  is  made  to  a  recent  text  book  by  the  PI  on  this  subject  [5]. 

4.1  Governing  Equations 

At  pressures  and  temperatures  normally  used  in  CVD  of  structural  materials,  the  mean 
fi’ee  path  of  the  gas  molecules  is  much  smaller  than  the  reactor  dimensions  and  as  such,  a 
continuum  analysis  is  adequate  to  study  the  transport  and  kinetics  inside  the  reactor.  The 
following  are  the  overall  mass,  momentum,  energy,  and  individual  species  mass  conservation 
equations,  respectively. 

^  +  Vv,u  =  0  (1) 

^  +  V-(pUU)  =  -VP+V-pl+pg  (2) 

^+V.(pTO)  =  -V.q+:|^+<I.  +  Q„  (3) 

+  f  (4) 

^  y=i 

where  p  is  total  density,  p,  is  density  of  species  / ,  U  is  mass-averaged  gas  velocity,  P  is  pressure, 

=> 

n  is  the  viscous  stress  tensor,  g  is  the  gravitational  vector,  h  is  enthalpy,  q  is  the  mean  heat  flux 
vector,  and  <I)  is  the  mean  flow  dissipation  rate.  Qext  includes  external  heat  sources  including 
radiation,  j;  is  the  diffusive  flux.  The  last  term  in  Eq.  (4)  represents  species  production  or  loss 
through  7  chemical  reactions. 

The  above  transport  equations  need  to  be  augmented  with  an  equation  of  state  such  as  the 
ideal  gas  law. 


PM 

RT 


(5) 


Here  M  is  the  mixture  molecular  weight  and  R  is  the  gas  constant.  Further  details  on  some  of  the 
key  terms  in  the  above  equations  are  given  below. 
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The  enthalpy  is  the  mixture  enthalpy  of  N  species. 


N 

h  =  (6) 

/=! 

where  (o,  is  the  mass  fraction  of  the  species  given  by  pjp.  The  molecular  heat  flux  q  is  the 

sum  of  the  convective  and  interdifiusional  contributions,  qc  and  qa  where 

q,  =  -kVT  (7) 

»=1 


where  xris  thermal  conductivity.  The  temperature  Tis  related  to  enthalpy  h  via 


^=Z<y, 


1=1 


(9) 


where  /»/,  is  heat  (enthalpy)  of  formation  of  species  /  at  temperature  Tp  ,  and  Cpi  is  the  specific 
heat  at  constant  pressure  for  species  /. 

4.2  Conjugate  Heat  Transfer 

For  problems  in  which  significant  heat  transfer  effects  are  present,  it  is  sometimes 
necessary  to  perform  a  coupled  analysis  in  which  the  fluid  flow  analysis  and  a  heat  conduction 
analysis  mutually  affect  one  another.  Generally  this  is  only  required  when  the  thermal  coupling  is 
strong  both  ways,  and  separate  analyses  would  fail  to  account  for  the  mutual  influence  of  the  solid 
medium  on  the  fluid  medium  and  vice  versa.  To  allow  for  a  full  range  of  choices  in  solving  heat 
transfer  problems,  a  conjugate  procedure  has  been  developed  and  applied  to  CVD  problems.  The 
choices  include  full  “real-time”  transient  coupling,  one-way  coupling,  and  full  “pseudo-time” 
coupling  aimed  at  reaching  a  steady-state  condition  as  rapidly  as  possible.  The  need  to  exchange 
data  between  disciplines  is  a  common  requirement,  whether  in  an  iterative  manner,  as  in  the  cases 
of  full  coupling,  or  as  in  a  one  time  event  for  the  case  of  one-way  coupling.  The  basis  of  the 
approach  is  the  sequential  coupling  of  an  arbitrary  number  of  related  fluid  and  solid  zones  to 
interact  with  one  another,  through  the  use  of  suitable  heat  flux  and  temperature  compatibility 
conditions.  The  schematic  below  illustrates  this  point. 
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While  heat  flux  and  temperature  are  calculated  in  both  the  fluid  and  the  solid,  only  one  of 
them  needs  to  be  passed  across  the  boundary  to  act  as  a  boundary  condition.  If  both  flux  and 
temperature  from  both  zones  were  exchanged  between  one  another,  it  would  be  impossible  to 
apply  both  as  boundary  conditions;  i.e.,  the  boundary  conditions  would  be  over-specified.  As  a 
result,  a  choice  must  be  made  regarding  which  medium  should  receive  the  flux  as  a  boundary 
condition  and  which  should  receive  temperature.  This  choice  is  usually  made  on  the  basis  of  the 
ratio  of  the  diffusive  time  scales  of  the  two  zones  and/or  their  thermal  conductivities.  These  two 
ratios  also  affect  the  solution  strategy  under  certain  circumstances.  The  issues  that  must  be 
discussed  are  as  follows: 

boundary  conditions 

solution  strategies  for  transient  and  steady-state  calculations 

Also  required  is  discussion  regarding  data  exchange  mechanisms  and  the  implementation 
of  the  program  in  a  distributed  computing  environment. 

Thermal  Boundary  Conditions:  While  the  conditions  applied  at  the  common  interface  between 
the  solid  and  fluid  zones  are  strictly  compatibility  conditions,  in  practice  these  are  applied  as 
boundary  conditions  in  each  of  the  analysis  codes.  In  that  sense  the  approach  is  no  different  fi’om 
the  approach  that  would  be  adopted  if  the  two  analyses  were  to  be  run  separately;  i.e.,  not  under 
the  automatic  control  of  a  computer  code.  When  performing  such  analyses,  it  is  normal  to  apply 
either  temperature  or  heat  flux  as  boundary  conditions,  but  not  both  on  the  same  boundary.  Of 
course,  combinations  of  flux  and  temperature  can  be  used,  but  they  must  be  applied  to  distinct 
surfaces. 
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For  convenience,  in  the  code  described  in  this  work  both  the  flux  and  temperature  were 
obtained  in  both  zones.  Having  both,  even  though  both  need  not  be  exchanged,  is  useful  for  the 
purpose  of  establishing  convergence  by  comparing  the  implied  flux  on  one  side  of  the  common 
boundary  with  the  flux  applied  on  the  other  side.  As  shown  schematically  earlier,  above,  the  flux 
is  calculated  using; 


-k„ 


= 


Tm  = 
where: 


is  the  local  arc  length  on  the  given  surface. 
m  the  subscript,  refers  to  any  medium, 
is  the  local  heat  flux  in  the  medium, 
is  the  thermal  conductivity  of  the  medium. 
7^  is  the  local  temperature  of  the  medium. 


From  a  finite  difference  point  of  view,  the  normal  derivative  appearing  in  the  flux  can  be 
evaluated  explicitly  using  any  one-sided  difference  formula  at  the  end  of  each  iterative  step  of  the 
calculation.  Both  temperature  and  flux  are  stored  as  functions  of  surface  arclength  to  enable 
interpolation  to  be  used  in  the  process  of  transferring  data  across  the  boundary  between  the  zones. 
This  is  an  important  issue,  because  in  the  most  general  cases  the  computational  grids  in  each  zone 
will  not  be  aligned  with  one  another  along  the  common  surface  and,  therefore,  it  is  necessary  to 
interpolate  to  allow  data  obtained  in  one  zone  to  be  applied  as  boundary  conditions  in  the  other.. 

Solution  Strategies  for  Transient  and  Steady-State  Operation;  Depending  upon  the  type  of 
calculation  being  performed,  the  strategy  for  solving  the  problem  may  change.  While  it  is  true 
that  every  problem  could  be  treated  as  though  it  were  a  true  physical  transient,  this  approach 
wastes  an  unacceptable  amount  of  time  if  a  steady  state  is  the  ultimate  objective.  It  therefore  pays 
to  investigate  the  possibility  of  accelerating  convergence  in  such  cases.  In  addition,  even  in  cases 
where  a  true  transient  is  desired,  one  medium  may  have  a  significantly  faster  transient  response 
than  the  other.  This  allows  for  real-time  in  the  slower  medium  and  “pseudo-time”  in  the  faster 
medium.  This  implies  that  since  the  overall  system  thermal  response  is  determined  by  the  slower 
medium  and,  relatively  speaking,  the  faster  medium  responds  instantly,  there  is  no  reason  to  run  a 
physical  transient  in  both.  The  aim  of  doing  so  is  to  reduce  the  computational  workload  and 
thereby  decrease  calculation  time  and,  ultimately,  design  cycle  time.  The  maximum  saving  would 
occur  if  the  fluid  were  to  have  the  faster  thermal  response,  since  methods  then  could  be  employed 
to  run  the  fluid  calculation  in  “pseudo-time”  and  converge  it  faster.  If,  however,  the  solid  was  an 
insulator  relative  to  the  fluid,  then  significant  savings  could  be  obtained. 
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4.3  Multicomponent  Diffusive  Flux  Model 

We  need  to  provide  a  description  of  the  term  j,  in  Eq.  (4).  This  consists  of  two 
contributions:  diffiision  due  to  concentration  gradient  (ordinary  diffusion)  and  due  to  temperature 
gradient  (Soret  or  thermal  diffusion).  The  ordinary  diffusive  flux  for  two  components,  /  and  j,  is 
given  by  Pick’s  law. 

h=-pDij^oii  (10) 

Here  is  binary  difflisivity  of  /  into  j.  If  the  carrier  gas  is  the  major  component  and  all  active 

species  are  small  fractions,  then  the  mixture  may  be  considered  dilute.  In  this  case,  each  active 
species  /  may  be  thought  of  as  diffusing  into  the  carrier  gas,  and  we  can  use  the  corresponding 
binary  diffusivity.  However,  when  a  carrier  gas  is  not  present  or  present  only  in  small  quantities 
or  happens  to  be  a  light  molecule  such  as  hydrogen  or  helium,  the  dilute  species  approximation  is 
not  appropriate;  we  must  then  consider  multicomponent  mass  transfer.  The  latter  is  described  by 
Stefan-Maxwell  equations;  however  in  practice,  an  effective  multicomponent  difihisivity  for 
species  /  into  the  mixture  is  defined,  and  expression  (10)  is  used  with  A  instead  of  Ay.  The 
definition  of  A  is: 


AT 


y=i  % 


(11) 


Here  x,  is  the  mole  fraction.  It  is  noted  that  diffusive  fluxes  of  all  species  must  sum  up  to  zero, 
N 

i.e.,  2  j,.  =  0;  only  then  all  species  conservation  equations  would  add  to  give  the  mass 
/  =  1 

continuity  Eq.  (1).  In  practice,  this  condition  is  realizable  with  the  form  of  Eq.  (10)  only  if  all 
diffiisivities  are  equal;  the  latter  is  rather  unrealistic.  While  different  approaches  may  be  used  to 
satisfy  the  above  constraint,  a  rigorous  procedure  yields 


j/ 


M 


^AV 


^  PjM^ 


Pl 

M 


N 


^PjM^ 


(12) 


as  shown,  for  example,  by  Ramshaw  [12,13].  Though  expression  (12)  is  significantly  more 
complex  than  in  expression  (10),  it  conserves  mass  and  yields  zero  net  diffusive  flux.  This 
approach  is  better  than  that  used  in  the  literature  and  in  many  commercial  codes,  where 
expression  (10)  is  normally  used.  Another  issue  of  concern  when  expression  (10)  is  used  instead 
of  (12)  is  that  usually  one  of  the  species  from  the  set  of  n  is  not  solved,  and  its  mass  or  mole 
fraction  is  obtained  by  subtracting  the  sum  of  the  other  n-1  species  from  unity.  This  is  not  only 
subject  to  round-off  errors,  but  also  may  result  in  an  “asymmetric”  situation,  wherein  the  final 
results  may  depend  on  which  species  is  the  odd  one  not  to  be  solved  for. 
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4.4.  Thermal  Diffusion  Model 


The  Soret  (thermal)  difEision  flux  is  given  by 

j(,r  =  -A,rV(ln7-)  (13) 

Here,  Dij  is  the  thermal  diffusion  coefficient  for  species  /.  Modelers  generally  rely  on  empirical 
relations  to  evaluate  Dij.  For  many  species  of  interest  to  us,  such  relations  are  not  available,  and 
hence,  we  incorporated  a  first-principles-based  procedure  in  our  model  and  code  for  this  purpose. 
This  procedure  is  based  upon  the  one  described  by  Ramshaw  [12,13],  Di^r  is  written  in  the  form 

%MjKj  (14) 

M  p  M  ^ 

where  M  is  the  molecular  weight  of  species  /,  and  where  AT,  is  a  thermal  diffusion  coefficient 
given  by  the  relation 


Here  Pg  is  approximated  by  the  expression 


(15) 


XjXfUjiTf 

P  =-P  •  J 


where 


and 


ri  = 


2kT 


_1_ 


rijcxy 


(16) 


(17) 


(18) 


In  these  relations,  /w,  is  the  mass  of  a  single  molecule  of  species  /,  w/  is  the  number  density 
of  species  /,  and 


mprij 

m,  +  ntj 


(19) 
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(20) 
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YiTj 

Yi+rj 


Cfi, 


at+aj 


] 


2 


(21) 


where  a,  is  the  molecular  collision  diameter  of  species  /  (and  Oy  is  the  total  cross-section  for  “/•/’ 
collisions).  Note  that  t,  represents  the  mean  time  between  collisions  of  molecules  of  species  /. 
The  form  of  Eq.  (16)  differs  from  that  given  in  Ref  12,  and  was  suggested  by  Ramshaw  [14], 
Finally,  note  that  the  form  of  Eq.  (14)  guarantees  that  the  sum  of  the  thermal  diffusion  fluxes  over 
all  species  is  zero. 

4.5.  Gas  Phase  Reactions  Model 

j 

The  last  term  ^  R^j  in  Eq.  (4)  represents  the  sum  of  the  rate  of  production  or 
/=' 

consumption  for  species  /  from  each  reaction  j.  In  thermal  reactions,  the  forward  reaction  rate  is 
given  by  the  usual  Arrhenius  form: 

Rate  =  kf[pQ}iY‘  where  kf  =  AT^  expj^-.£  /  TJrj  (22) 

Here  p  is  the  mass  density  identified  in  Eq.  (1).  is  the  forward  reaction  rate  constant  which  has 
an  exponential  dependence  on  temperature;  E  is  the  activation  energy;  A  is  an  Arrhenius  pre¬ 
exponential  factor;  is  the  order  of  reaction  with  respect  to  species  /.  For  reversible  reactions, 
the  reverse  reaction  rate  constant  is  related  to  the  equilibrium  constant  by  k^  =  k^jK^  where  ATc  is 

the  equilibrium  rate  constant.  The  latter  can  be  computed  from  thermodynamic  properties  (Gibb’s 
free  energy)  using  standard  procedures.  We  use  notations  to  describe  and  procedures  to  handle 
chemical  reactions  similar  to  those  in  the  CHEMKIN  procedure  [15]. 


4.6  Surface  Process  Model 


Surface  reactions  involve  both  gas  phase  species  and  surface  products.  The  notation  and 
handling  procedure  for  surface  species  and  reactions  can  be  done  in  a  manner  similar  to  the  gas 
phase  reactions.  This  has  been  done  in  the  Phase  I  work.  Surface  activity  consists  of  a  sequence 
of  elementary  steps;  (1)  diffusion  of  reacting  species  to  the  surface,  (2)  adsorption  of  these  onto 
the  surface,  (3)  surface  chemical  reaction,  and  (4)  desorption  of  the  products.  The  proper 
approach  is  to  adequately  compress  the  available  information  on  any  or  all  of  the  above  four  steps 
(in  practice,  usually  one  of  them  would  be  the  rate  controlling  step  in  many  situations)  and  the 
detailed  surface  chemistry  into  appropriate  boundary  conditions  [5].  Computation  of  forward, 
reverse  and  equilibrium  surface  rate  constants  is  done  in  a  manner  identical  to  that  for  a  gas  phase 
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reaction.  One  major  difference  is  the  need  to  include  surface  site  as  an  active  participant  in  a 
reaction  so  that  the  rate  would  depend  on  the  availability  of  active  sites. 

4.7  Transport  and  Thermodynamic  Properties 

For  the  deposition  of  a  given  material,  the  transport  properties  such  as  viscosity,  thermal 
conductivity  and  diffusivity  must  be  specified  as  a  function  of  temperature  and  pressure.  These 
are  obtained  fi'om  Lennard-Jones  parameters  for  each  of  the  chemical  species.  Note  that  all 
physical  properties  vary  across  the  reactor  with  a  variation  in  temperature  and  composition; 
therefore  assumption  of  constant  properties  is  not  appropriate.  Thermodynamic  properties  such 
as  heat  capacity,  enthalpy,  entropy,  and  fi-ee  energy  for  each  species  is  obtained  fi'om  the  NASA 
Lewis  database.  In  the  current  work,  we  have  obtained  a  Sandia  data  base  and  added  this  data 
base  to  our  code.  This  database  contains  data  for  many  silane  and  chlorosilane  species  used  in 
structural  materials  preparation. 

4.8  Boundary  Conditions 

A  complete  problem  definition  requires  specification  of  boundary  conditions.  At  a 
subsonic  inflow  boundary,  total  pressure,  total  temperature,  and  flow  angles  are  specified.  This 
has  proven  very  effective  for  a  wide  variety  of  cases.  For  species,  Danckwert's  boundary 
condition  is  used  at  the  inlet,  which  equates  the  incoming  mass  flux  to  the  net  convective  and 
diffusive  flux  at  the  inlet.  This  is  an  alternative  to  specifying  the  mass  fractions  at  the  inlet. 
Outflow  boundary  conditions  require  setting  downstream  pressure  and  extrapolation  of  species 
densities  and  velocity. 

At  the  walls,  a  no-slip  condition  is  imposed  for  velocities.  At  the  substrate  and  boundaries 
where  deposition  occurs,  a  mass  balance  for  each  species  equates  the  flux  to  the  growth  surface 
and  the  net  production  of  that  species  in  surface  reactions.  In  representing  the  surface  reactions, 
one  needs  to  include  detailed  mechanisms  including  adsorption  and  desorption  kinetics,  as 
discussed  above  [5].  Note  also  that  the  net  growth  rate  results  in  a  non-zero  mass-averaged 
velocity  at  the  growth  surface. 

For  the  energy  equation  either  the  temperature  or  heat  flux  at  the  reactor  wall  can  be 
specified.  It  is  well  known  that  in  many  reactors  the  wall  and  substrate  temperatures  are  not 
constant;  nevertheless,  the  constant  temperature  boundary  condition  has  been  routinely  used  in 
the  literature.  Realistically,  the  heat  transfer  within  the  wall  (or  wafer)  must  be  considered  with 
proper  matching  to  the  heat  flow  within  the  reactor  on  one  side  and  the  surroundings  on  the  other 
side.  The  time  evolution  of  temperature  within  a  solid  for  constant  properties  is  given  by 

=  (23) 

Here  the  subscript  s  stands  for  solid.  Q  denotes  heat  generation  within  the  solid  (due  to  inductive 
heating,  for  example).  The  boundary  conditions  on  the  side  of  the  reactor  volume  must  include 
heat  flow  due  to  both  convection  and  radiation.  The  radiative  heat  transfer  term  is  necessary 
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because  of  the  radiative  heat  exchange  between  heated  susceptor,  walls  and  windows.  For  the 
outer  wall,a  convective  description  may  be  adequate.  The  above  analysis,  called  conjugate  heat 
transfer,  is  a  must  for  high  temperature  reactor  operation  to  accurately  describe  wall  and 
susceptor  temperature  distributions  and  has  been  described  in  Section  4.2.  Note  that  in  the  case 
of  deposition  on  fibers,  the  substrate  fiber  moves  continuously  and  heat  transfer  within  this 
moving  boundary  must  be  considered. 

4.9  Solution  and  Numerical  Procedure 


The  basic  numerical  procedure  and  in-house  computational  fluid  dynamics  code  has  been 
developed  by  SRA  personnel  over  many  years  through  core  funding  from  Air  Force  and  NASA 
agencies.  Our  code  has  been  applied  to  a  variety  of  fluid  flow  problems  by  both  us  and  personnel 
from  the  above  two  agencies.  The  current  capability  includes  two  and  three  dimensional  analysis, 
reacting  flow,  two  phase  flow,  and  multiple  turbulence  models.  We  also  have  grid  generation 
capability  to  analyze  complex  3-D  nonorthogonal  geometries.  This  capability  is  essential  in  the 
proposed  work  which  will  involve  modeling  of  3-D  objects.  Applications  of  our  code,  numerical 
procedure  and  grid  capability  have  been  made  to  a  wide  range  of  practical  problems  for  DoD  and 
NASA  including  turbomachinery,  primary  and  secondary  gas  paths,  re-entry  vehicles,  rocket 
motors,  inlets,  nozzles,  seals,  chemical  lasers,  combustion,  two  and  multi-phase  flows,  crystal 
growth,  CVD,  and  reactive  ion  etching.  This  code  was  adapted  in  Phase  I  with  specific 
modifications  described  earlier  for  CVD  reactor  analysis. 

The  numerical  algorithm  which  forms  the  basis  for  this  code  is  described  below.  The 
procedure  to  solve  the  governing  equations  is  a  consistently  split  linearized  block  implicit  scheme 
originally  developed  by  Briley  and  McDonald  [16]  at  SRA  and  embodied  in  a  computer  code 
termed  MINT.  The  basic  algorithm  has  been  further  developed  and  applied  to  both  laminar  and 
turbulent  flows  (see  Briley  and  McDonald  [17]).  The  method  can  be  outlined  as  follows:  the 
governing  equations  in  a  dimensionless  form  are  replaced  by  an  implicit  time  difference 
approximation.  Terms  involving  nonlinearities  at  the  unknown  time  level  are  linearized  by  Taylor 
series  expansion  about  the  solution  at  the  previous  known  time  level,  and  spatial  difference 
approximations  are  introduced.  The  result  is  a  system  of  multi-dimensional  coupled  (but  linear) 
difference  equations  for  the  dependent  variables  at  the  unknown  or  implicit  time  level.  To  solve 
these  difference  equations,  the  Douglas-Gunn  procedure  for  generating  alternating  direction 
implicit  (ADI)  splitting  schemes  is  used.  This  ADI  splitting  technique  leads  to  systems  of  coupled 
linear  difference  equations  having  narrow  block-banded  matrix  structures  which  can  be  solved 
efficiently  by  standard  block-elimination  methods. 

Two  critical  elements  of  a  numerical  algorithm  are:  (a)  accuracy  and  (b)  stability  and 
convergence  rate.  Regarding  accuracy,  most  of  the  finite-difference  algorithms  developed  for 
complex  fluid  dynamic  problems  utilize  either  an  upwind  scheme  or  a  central  difference  scheme 
with  artificial  dissipation  terms.  Artificial  dissipation  terms  in  a  central  difference  scheme  serve  to 
provide  a  degree  of  upwinding  to  a  central  difference  scheme,  depending  upon  the  magnitude  of 
the  terms.  In  turn,  the  artificial  dissipation  terms  reduce  accuracy  of  the  central  difference 
scheme.  The  use  of  an  upwind  scheme  or  a  central  difference  scheme  with  artificial  dissipation 
has  advantages  and  disadvantages  and  is  a  subject  of  active  research.  The  advantage  of  a  central 
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difference  scheme  with  artificial  dissipation  is  that  the  magnitude  of  the  artificial  dissipation  can  be 
controlled.  For  example,  in  viscous  regions  of  the  flow  or  for  low  Reynolds  number  flows  the 
dissipation  terms  can  be  significantly  reduced  or  eliminated,  providing  better  accuracy  than  an 
upwind  scheme  on  the  same  finite  difference  stencil.  This  is  the  strategy  adopted  in  MINT.  The 
mathematical  form  of  the  dissipation  terms  utilized  is  based  on  the  convective-diffusive  character 
of  the  Navier-Stokes  equations  for  single  species  transport.  While  this  form  of  the  terms  is 
suitable  for  the  CVD  problem,  they  may  not  be  optimal.  For  the  CVD  problem,  species 
generation  through  chemical  reactions,  heat  transfer,  and  buoyancy  significantly  affect  flow 
evolution  and  these  effects  should  be  included  in  choosing  a  form  for  the  artificial  dissipation 
terms.  This  aspect  has  been  considered  in  the  present  work. 

Regarding  convergence  rate  of  the  numerical  algorithm,  matrix  conditioning  techniques 
have  been  utilized  in  MINT  to  accelerate  convergence  to  steady  state  when  only  the  steady  state 
solution  is  of  interest.  These  matrix  conditioning  techniques  spatially  alter  the  time-step  based  on 
the  local  eigenvalues  of  the  coefficient  matrices  of  the  convective  and  diffusive  terms.  Again,  for 
the  CVD  problem,  there  are  additional  processes  of  generation  through  chemical  reactions,  heat 
transfer,  and  buoyancy  that  need  to  be  included  in  the  matrix  conditioning  procedures.  Another 
element  that  affects  the  numerical  algorithm  is  that  for  very  low  Mach  number  flows  such  as  those 
occurring  in  CVD-related  problems,  the  numerical  solution  of  the  equations  (l)-(3)  with  the 
equation  of  state  P/p  -  RT  leads  to  two  problems,  viz.  (i)  slow  convergence,  and  (ii)  reduced 
time  accuracy,  due  to  the  ill-conditioned  nature  of  some  of  the  coefficient  matrices.  These 
problems  can  be  overcome  by  replacing  the  equation  of  state  by  P^lp  =  RT,  where  Po  is  a 
constant  pressure  (which  corresponds  to  the  almost  constant  pressure  in  the  system  when  the 
Mach  number  is  very  low),  and  by  solving  the  system  of  equations  (l)-(3)  for  the  dependent 
variables,  P,  U,  and  h. 

The  solution  efficiency  of  species  equations  depends  on  the  time  scale  for  the  chemical 
reactions.  For  example,  in  our  previous  study  on  projectile  base  region  combustion  flows  [18], 
the  reaction  terms  were  found  to  yield  a  stiff  system  of  equations  due  to  widely  varying  time 
scales  in  the  problem.  The  implicit  scheme  discussed  above  was  successfully  applied  by  fully 
coupling  key  species  equations  with  the  flow  and  energy  equations.  We  also  developed  a  matrix 
conditioning  technique  to  allow  treatment  of  widely  varying  (time-scale)  source  terms.  This 
technique  employs  an  ad  hoc  spatially  varying  conditioning  factor  in  the  species  equations  which 
is  chosen  to  insure  diagonal  dominance  of  the  coupled  block  matrices  at  each  grid  point.  The 
approach  is  equivalent  to  rescaling  the  time  step  operators  for  the  density  (p)  and  species  mass 
fractions  (cy,)  in  each  of  the  coupled  species  equations.  This  is  a  justifiable  operation  as  long  as 
the  resulting  system  of  difference  equations  converges  satisfactorily,  as  was  demonstrated  by  Ref. 
18.  For  transient  problems  the  iterative  procedure  discussed  below  permits  the  use  of  scaling 
techniques  during  the  iteration  at  each  time  step. 

The  solution  procedure  described  above  can  be  used  to  obtain  both  steady-state  and  time- 
accurate  unsteady  flow  solutions.  The  convergence  to  a  steady-state  solution  can  be  sped  up  by 
using  matrix  preconditioning  techniques,  which  effectively  lead  to  the  use  of  different  “time”  steps 
for  each  grid  point.  Unsteady  flow  solutions  can  be  obtained  in  two  ways:  (a)  by  choosing  a 
sufficiently  small  time  step  (without  using  matrix  preconditioning),  or  (b)  by  using  an  iterative 
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procedure  to  obtain  the  solution  at  the  next  time  step.  This  latter  approach  allows  the  use  of 
larger  time  steps  without  reducing  the  time  accuracy,  especially  if  the  predominant  errors  are  due 
to  the  ADI  factorization,  at  the  cost  of  a  few  iterations  per  time  step.  This  may  be  advantageous 
for  problems  with  multiple  time  scales. 


4.10  Lower  Order  Models 


As  mentioned  earlier,  simulations  of  materials  processing  problems  typically  involve  a 
large  number  of  chemical  species  in  a  large  number  of  reactions.  As  each  of  these  species  requires 
a  conservation  equation  in  the  equation  set,  a  typical  problem  consists  of  seeking  multidimensional 
solutions  to  a  large  number  of  equations.  Given  the  stiffness  of  the  reaction  terms,  it  is  not 
unusual  to  experience  “slow  convergence”;  this,  combined  with  the  size  of  the  matrix  due  to  a 
large  number  of  equations  encountered  in  materials  processing  simulations,  makes  turnaround 
times  of  a  few  days  a  common  experience.  While  computational  procedures  to  accelerate 
convergence,  combating  stiff  reaction  terms  through  proper  time  scaling,  etc.,  command  the 
attention  of  researchers,  and  we  also  focus  on  such  procedures  in  the  Phase  I  and  an  anticipated 
Phase  II,  it  is  also  useful  to  examine,  a  priori,  the  species  and  reaction  set  carefully.  It  is  possible 
to  identify  trace  species  and  reactions  too  slow  in  the  time  scale  of  interest.  While  intuition  can 
guide  in  this  identification  task,  it  may  be  more  complicated  in  a  flow  environment.  Keeping  this 
in  mind,  one  can  use  “lower  order”  models  to  achieve  this  goal.  By  this  we  mean  0  or  1 
dimensional  models,  in  contrast  to  the  2  and  3-d  models. 

A  1-d  model  consists  of  radially-averaging  the  governing  equations  presented  earlier  and 
obtaining  transport  equations  in  the  flow  direction.  The  activities  in  the  radial  wall  are  self- 
consistently  converted  to  source  terms  in  the  resultant  equations  through  the  application  of  radial 
boundary  conditions.  A  0-d  model  consists  of  volume-averaging  the  governing  equations 
presented  earlier  and  obtaining  algebraic  conservation  equations  for  mass  and  energy.  Reaction 
terms  in  the  gas  phase  are  intact,  while  the  surface  reaction  terms  are  folded  properly  into  the 
analysis  by  weighting  them  with  the  surface-to-volume  ratio. 

While  these  lower  order  models  do  not  capture  the  entire  picture  of  a  complex  process  in 
a  complex  geometry,  the  idea  of  using  them  is  not  as  an  alternative  to  multidimensional 
simulations,  but  rather  as  a  preliminary  or  screening  procedure.  Though  such  a  procedure  is  not 
common  in  traditional  CFD  of  subsonic  and  supersonic  analyses  and  turbomachineiy  studies,  it 
makes  a  lot  of  sense  in  chemically  reacting  flows  common  in  materials  processing.  In  the  present 
work,  existing  SRA  0-d  and  1-d  codes  have  been  adapted  to  apply  to  CVD  problems.  The 
analysis  has  been  repeatedly  used  in  an  effort  to  come  up  with  a  “manageable  chemistry  set”  that 
contains  all  aspects  of  the  processing. 
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5.  Preform  Analysis:  Chemical  Vapor  Infiltration  Model 

5.1  Modeling  Approach 

CVI  analysis,  in  general,  has  two  components: 

(1)  transport  and  reactions  external  to  the  preform  in  the  reactor  and 

(2)  transport  and  kinetics  inside  the  preform. 

These  are  coupled  and,  perhaps,  can  be  handled  in  a  “conjugate  mass  transfer”  approach 
analogous  to  conjugate  heat  transfer  discussed  earlier.  That  is,  the  coupling  occurs  at  the 
boundaries.  Note  then  that  item  (1)  regarding  the  transport  external  to  the  preform  is  similar  to 
CVD  analysis.  Item  (2)  requires  model  development,  which  is  the  subject  of  this  chapter. 
Coupling  issues  will  be  addressed  in  Phase  II. 

In  principle,  the  transport  inside  the  preform  is  governed  by  the  same  equations  described 
in  the  previous  section;  however,  there  are  several  important  points  worth  noting  here:  (1)  In 
isothermal  no-flow  preforms,  all  the  terms  corresponding  to  flow  and  flow  velocity  U  cm  be 
dropped.  (2)  In  forced  flow  reactors,  the  flow  may  be  represented  by  Darcy’s  law.  (3)  Terms 
such  as  dissipation,  O  in  Eq.  3,  are  not  important.  (4)  The  term  Qext  in  Eq.  3  represents  the 
heating  mode  (rf  volume  heating,  induction  heating,  microwave  heating,  etc.)  and  each  mode  must 
be  appropriately  modeled  as  source  terms.  If  microwave  heating  is  involved.  Maxwell’s  equations 
need  to  be  solved.  (5)  The  reaction  terms  Ry  in  Eq.  4  corresponds  to  gas  phase  reactions  inside 
the  pore.  The  species  /  may  be  a  parent  gas,  intermediate,  precursor  or  material  to  be  deposited. 
Note  that  gas  phase  reactions  inside  the  pores  may  not  be  important  but  it  is  hard  to  make  such  a 
generalization;  it  is  strongly  system  dependent  and  so  we  include  this  term.  (6)  Compared  to  the 
time  of  pore  evolution,  the  time  scales  for  mass,  momentum  and  energy  transfer  are  fast  and 
therefore  steady  state  forms  of  the  above  equations  are  adequate. 

5.2  Pore  Diffusivitv 

Next  let  us  discuss  the  nature  of  diffusion  inside  the  preform.  The  pore  diffosivity  in 
general  may  not  be  governed  by  molecular  diffusion.  The  diffiision  mode  is  determined  by  the 
pressure  and  pore  radius.  When  the  pore  radius  becomes  comparable  to  the  mean  free  path  X, 
then  molecules  collide  with  pore  walls  and  rebound;  collision  with  other  molecules  is  insignificant 
compared  to  collision  against  pore  walls.  Diffusion  under  these  conditions  is  called  Knudsen 
diffusion  and  given  by 


=  4850c/ 
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for  A  >  \Qd 


(24) 


Here  d  is  the  pore  diameter.  Mi  is  the  molecular  weight  of  species  /,  and  Tis  gas  temperature. 
For  0.01c/  >X,  ordinary  diffusion  applies  and  diffusivities  can  be  calculated  using  Lermard-Jones 
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parameters  and  well-known  formulae.  In  the  transition  region,  an  effective  difiusivity  comprising 
of  the  two  methods  can  be  used. 

5.3  Pore  Evolution 


The  governing  transport  equations  when  solved  can  provide  local  species  density, 
temperature,  etc.  Avithin  the  preform.  The  time  evolution  of  the  accessible  porosity  and  its  spatial 
dependence  within  the  preform  (i.e.,  microscopic  pore  evolution)  can  be  obtained  from 
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Here  s  is  the  porosity,  Md  and  pd  are  molecular  weight  and  density  of  deposit,  Rk  is  the  of « 
surface  and/or  gas  phase  reactions  which  produces  or  consumes  the  deposit  material,  is  the 
corresponding  stoichiometric  coefficient.  This  dynamic  equation  for  porosity  evolution  must  be 
solved  along  with  the  governing  equations. 

5.4  Network  Representation 

This  brings  us  to  the  final  question  with  respect  to  preform  modeling.  For  a  simulation  to 
be  useful,  we  must  take  into  account  the  preform  architecture  in  a  meaningful  way.  So,  how  do 
we  model  the  porous  network?  Or,  stated  differently,  what  geometrical  coordinates  do  we  use  to 
write  and  solve  our  equations? 

One  may  assume  that  the  preform  consists  of  uniform  and  unidirectional  cylindrical  pores. 
This  is  an  extremely  idealistic  but  simplistic  approach.  This  has  been  attempted  in  the  literature 
[19]  and  it  is  not  known  if  the  results  reproduce  observed  experimental  behavior  (no  comparison 
was  made).  We  do  not  believe  that  this  approach  would  provide  a  meaningful  representation  of 
any  of  the  preforms  used  in  CVI.  There  have  been  attempts  to  represent  the  random  porous 
network  using  different  coordinate  representations  [9,1 1].  Such  representations  may  well  be 
preform-specific  or  material  specific  and  no  detailed  comparison  with  experiments  is  available  to 
assess  the  validity  of  these  representations.  An  elegant  approach  to  model  transport  through  a 
network  of  fibers  inside  the  preform  is  called  the  ‘dusty  gas  model’  and  used  successfully  in  the 
field  of  transport  phenomena  in  porous  catalysts.  Note  that  the  gas  flow  through  randomly 
packed  catalyst  particles  is  similar  to  the  flow  of  gas  phase  precursors  through  the  porous 
network.  The  heterogeneous  catalysis  area  in  chemical  engineering  literature  has  abundant 
examples  of  this  approach  and  a  nice  discussion  can  be  found  in  a  text  book  by  Jackson  [20].  We 
recently  introduced  this  technique  in  the  field  of  CVI  and  has  further  developed  it  under  the 
present  program.  The  preform  architecture  consists  of  randomly  oriented  fibers  in  a  network.  In 
this  approach,  the  flow  through  a  network  of  porous  media  is  thought  to  be  equivalent  to  the  flow 
of  gas  on  a  molecular  scale  through  an  assembly  of  stationary  obstacles  dispersed  in  the  gas 
(hence  the  name  ‘dusty-gas’  model).  Based  on  this  premise,  there  has  been  a  systematic 
development  of  multicomponent  species  flux  relations  for  transport  within  the  porous  network; 
the  development  includes  pore  diffusion  (both  molecular  and  Knudsen  diffusion)  and  viscous 
flow.  The  dusty  gas  model  gives  expressions  for  molar  flux  of  each  species,  in  a  random  porous 
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network,  as  a  function  of  mole  fractions,  diffusivities,  transport  parameters,  and  gradients  in 
composition  and  pressure.  The  resultant  expressions  for  flux  in  a  multicomponent  mixture  are  too 
wieldy  to  reproduce  here  and  are  as  given  in  Ref  [20],  These  expressions  need  to  be  used  along 
with  the  governing  equations  presented  earlier. 

5.5  Preform  Heating 

The  preform  can  be  heated  by  any  one  of  several  techniques.  Large  pieces  are  heated  in 
high  temperature  furnaces.  We  have  previously  shown  [21]  that  microwave  heating  can  be  used 
to  heat  small  pieces,  whereby  it  is  possible  to  get  the  “inside-out”  heating  and  densification.  We 
have  also  investigated  pulsed  volume  heating  techniques  [22].  It  is  our  objective  to  investigate 
various  heating  techniques  in  Phase  I  and  II,  in  order  to  obtain  the  desirable  infiltration 
characteristics.  Indeed,  it  is  one  of  the  attractions  of  a  model  where,  cost-effectively  and  quickly, 
one  can  investigate  several  heating  techniques  and  designs  prior  to  the  effort  to  build  and  test 
them.  Conventional  heating  techniques,  from  a  modeling  point  of  view,  require  only  specification 
of  heat  flux  and/or  temperatures  at  the  boundaries.  Using  these  conditions,  transport  equations 
are  solved  to  obtain  temperature  profiles  within  the  preform.  On  the  other  hand,  if  microwave 
heating  is  used,  the  situation  is  more  complex,  and  we  addressed  this  issue  in  Phase  I.  As  a  result, 
we  added  another  module  to  our  model  and  code  that  can  investigate  microwave  heating. 

Heat  dissipation  in  the  ceramic  depends  on  the  magnitude  of  the  local  electric  field.  This 
is  found  from  Maxwell’s  equations,  which  for  the  system  of  interest  reduce  into  a  single  phasor 
expression,  as  follows  [23]: 

V^E,-r^E,=0  (26) 

in  which  Es  is  the  electric  field.  Quantity  y  is  known  as  the  complex  propagation  constant  and  is 
given  by 


y  =  o)  -yJjUQSQ  {-k'+ik")  =  a  + i/3 


(27) 


where  the  attenuation  and  phase  factors  are  given  by  the  following  expressions,  respectively: 
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In  the  above  equations,  co  (=2n{)  is  the  angular  frequency  of  the  applied  field, 

/c’  (=  sVso  =  oyco '  ^/soj  is  the  relative  loss  factor.  The  relative  dielectric  constant  accounts  for  the 
polarization  of  a  medium  caused  by  an  applied  electric  field.  The  relative  loss  factor  accounts  for 
the  development  of  currents  that  lead  to  the  thermal  dissipation  of  energy. 

To  determine  the  power  dissipated  in  the  composite,  O,  we  consider  the  propagation  of 
plane  waves  in  the  z-direction,  incident  symmetrically  onto  the  sides  of  the  composite  located  at 
z  =  ±L.  Invoking  appropriate  continuity  conditions  at  the  interfaces  between  ceramic  and 
surrounding  medium,  one  can  obtain  an  expression  for  Ex 


,-L<z<L 


(30) 


where  subscripts  zero  and  c  denote  ffee-space  and  ceramic  properties,  respectively.  The  intrinsic 
impedance  of  medium  m  can  be  obtained  from 


(31) 


The  power  dissipated  per  unit  volume  of  composite  can  then  be  obtained  from 


<t>  =  ^asak'El 


(32) 


In  the  last  equation,  the  relative  loss  factor  A”  is  a  function  of  porosity.  As  the  ceramic  absorbs 
energy  and  its  temperature  increases,  the  loss  factor  E'  increases  as  well,  and  even  more  energy  is 
absorbed.  This  feedback  mechanism  can  lead  to  thermal  runaway  [24].  The  simplest  form  of 
porosity  dependence  is 

k  ^{\-(l))kl+(l)kl  (33) 

where  ^  is  the  total  porosity. 
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5.6  Thermal  Stresses 


During  CVI  the  ceramic  composite  is  subjected  to  high  temperatures  that  may  lead  to  the 
development  of  significant  thermal  stresses.  In  general,  thermal  stresses  may  be  induced  by  a 
nonuniform  temperature  variation  within  the  medium,  differences  in  thermomechanical  properties 
in  a  multicomponent  medium,  or  by  directional  variations  in  the  properties  of  the  material.  In  this 
study  our  purpose  is  to  assess  the  thermal  stresses  due  to  temperature  gradients  to  assure  that 
crack  initiation  does  not  occur. 

A  simple  technique  is  presented  below  and  a  more  detailed  study  can  be  included  in  Phase 
n.  To  evaluate  the  stress  distribution  we  consider  that  the  thickness  of  the  slab  is  small  compared 
to  the  lateral  dimensions,  the  slab  faces  are  free  of  tractions,  and  that  the  system  is  symmetric 
about  the  midplane.  This  leads  to  a  plane  stress  problem  [25],  i.e.. 


^  zz  ~  ^zy  ~  ^zx  ~  ^  xy  ~  ^ 


(34) 


The  stresses  can  then  be  readily  found  from  the  compatibility  conditions 

1  ri 


UyE 

<^xx=(^yy=(^  =  7^ 

1-  y 


jja  T(t,z)dz-T(t.z) 


(35) 


whereas  the  strains  are  obtained  from  the  stress-strain  relations 


^xx  -  ^yy  - 


(36) 


a-! 


\-v 


~ T z)dz  +  {\+v)T{t,z)-{\-v) To 


(37) 


'zy 


'ZX 


"^xy  ~  (^ 


(38) 


The  thermomechanical  properties  are  the  coefficient  of  thermal  expansion,  ay.  Young’s  modulus 
of  elasticity,  E,  and  Poisson’s  ratio,  v.  These  properties  are  assumed  to  be  constant.  It  should  be 
noted  that  the  thermal  stress  problem  can  be  solved  independently  from  the  temperature 
distribution.  This  is  because  the  thermoelastic  dissipation  term  in  the  energy  balance  is  negligible. 
In  Eqs.  (36)  and  (37),  To  is  a  base  temperature  at  which  the  strains  in  the  material  are  zero  (here 
taken  as  300  K).  We  realize  that  the  above  equations  apply  to  a  homogeneous  isotropic  medium. 
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However,  we  still  use  these  equations  to  obtain  an  order  of  magnitude  estimate  of  the  stresses  in 
the  composite. 

5.7  Numerical  Method 


The  model  consists  of  a  set  of  coupled  partial  differential  equations  for  the  gaseous  species 
mole  fractions,  species  fluxes,  temperature,  and  porosity  subject  to  the  appropriate  boundary  and 
initial  conditions.  The  direct  solution  of  this  set  of  equations  using  the  method  of  lines  would 
result  in  a  mixed  differential-algebraic  system.  Alternatively,  one  can  obtain  an  equivalent  system 
consisting  only  of  differential  equations.  This  reduction  in  index  can  be  obtained  by  inverting 
numerically  at  every  time  step  the  system  of  linear  equations  represented  by  the  Dusty  Gas 
multicomponent  flux  equations.  This  procedure  reduces  the  number  of  dependent  variables  and 
increases  the  efficiency  of  the  method. 

The  discretization  of  the  spatial  derivatives  was  accomplished  using  orthogonal 
collocation  on  finite  elements  with  B-splines  basis  functions  [26].  The  number  of  collocation 
points  was  chosen  such  that  the  spatial  and  temporal  variations  of  the  power  dissipated  within  the 
composite  could  be  captured.  The  resulting  set  of  ordinary  differential  equations  was  integrated 
in  time  using  a  variable-step  variable-formula  method.  This  method  has  been  used  in  all  our 
previous  CVI  investigations  [21,  22,  27]. 
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6.  CVD  Results 


In  this  section,  we  discuss  results  from  several  demonstration  computations  of  CVD. 

First,  computations  specifically  performed  to  validate  the  model  are  discussed  in  Section  6.1.  We 
chose  low  pressure  CVD  of  silicon  for  this  purpose,  simply  because  it  is  an  extremely  well- 
characterized  and  understood  system.  In  addition,  these  calculations  demonstrate  the  generality 
and  application  of  our  models  and  code  for  electronics  materials,,  which  are  also  of  interest  to  the 
Air  Force  and  command  a  significant  market  for  the  software  from  Phase  III  efforts.  Next,  in 
Section  6.2,  we  discuss  CVD  of  high  temperature  silicon  carbide.  Two  geometries  are 
considered:  a  tubular  reactor  and  a  conventional  deposition  reactor  with  a  flat  substrate.  Section 
6.3  focuses  on  boron  deposition  on  tungsten  wires.  It  is  emphasized  that  the  Phase  I  scope  is 
limited  and  other  materials  of  interest,  for  example,  other  carbides  and  nitrides,  will  be  considered 
in  a  Phase  II  effort.  For  a  general  outline  of  Phase  II  directions,  please  see  Section  8. 

6.1  Model  Validation;  Silicon  CVD 


Our  major  objective  in  performing  these  computations  is  to  validate  the  thermal  diffusion 
model  and  computation  of  thermal  diffusion  coefficients  using  the  Ramshaw  formulation 
discussed  in  Section  4.4.  Kleijn,  et  al.  [28]  conducted  a  detailed  study  of  low  pressure  CVD  of 
silicon  wherein  they  used  empirical  coefficients  for  thermal  diffusion.  They  considered  two 
systems:  (a)  silane  in  a  carrier  gas  of  H2;  molecular  weights  of  silane  and  H2  are  about  15:1, 
(b)silane  and  H2  in  a  carrier  gas  of  N2;  silane  and  N2  molecular  weights  are  comparable.  Note  that 
thermal  diffusion  depends  on,  in  addition  to  existing  temperature  gradient,  the  molecular  size  and 
weight  of  the  components.  Here  we  repeated  Kleijn’s  computations  using  his  empirical 
coefficients  and  also  the  Ramshaw  model  from  Section  4.4.  A  detailed  comparison  is  made  to 
validate  the  thermal  diffusion  model  and  in  general,  the  CVD  model  and  our  code.  This  is 
important,  since  many  of  the  chemistries  used  for  emerging  high  temperature  materials  do  not 
have  empirical  data  for  transport  properties. 

The  schematic  of  the  CVD  reactor  from  ref  [28]  is  shown  in  Fig.  1 .  Gases  enter  through 
an  inlet  section  at  the  top,  come  down  the  metal  pipe,  flow  through  the  weir  and  over  the 
substrate,  and  down  through  an  exhaust  pipe.  The  geometry,  as  shown,  is  cylindrically  symmetric 
and  hence,  the  computations  are  2-d  axisymmetric.  The  outer  walls  are  at  room  temperature.  We 
treated  the  metal  pipe  in  two  different  ways:  (a)  wall  temperature  fixed  at  290  K,  and  (b)  heat 
transfer  through  the  wall,  using  the  conjugate  analysis  described  in  Section  4.2.  The  wafer 
temperature  is  taken  to  be  1,000  K,  as  in  ref  [28].  Though  the  heat  transfer  through  the  wafer 
can  be,  in  principle,  analyzed  using  the  conjugate  model,  it  was  not  done  since  ref  [28]  does  not 
provide  a  heat  flux.  The  area  surrounding  the  wafer  and  the  inside  wall  of  the  exhaust  pipe  are 
taken  to  be  adiabatic.  Figure  2  shows  the  grid  used  in  this  study.  It  consists  of  53  grid  cells  in  the 
axial  direction  and  38  in  the  radial  direction.  Critical  regions,  such  as  the  substrate  and  wall 
regions,  have  dense  packing  of  grid  cells  relative  to  other  places. 

For  our  purpose  here,  a  simple  chemistry  model  is  used  for  silicon  deposition  from  silane. 
As  in  ref  [28],  gas  phase  reactions  are  assumed  negligible  at  low  pressures.  (Note  that  we 
demonstrate  our  capability  to  analyze  transport  of  a  large  number  of  species  and  chemical 
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reactions  in  the  next  section  for  a  problem  of  relevance  to  this  Phase  I  program;  namely,  SiC 
deposition  for  aerospace  applications.)  In  this  section,  our  objective  is  to  compare  against  results 
from  ref  [28]  and  validate  our  code.  The  deposition  is  due  to  surface  adsorption  of  silane, 
followed  by  the  decomposition  of  silane  into  silicon  and  hydrogen,  leading  to  an  overall  reaction; 

SiH4(^)^Si(s)+2H,te) 

Kleijn,  et  al.  [28]  postulate  that  surface  adsorption  of  silane  is  the  rate-limiting  mechanism,  in 
which  case,  the  deposition  rate  can  be  given  by: 

The  values  for  the  constants  above  and  reaction  conditions  are  given  in  Table  I 

Table  I:  Si  CVD  Simulation 


Pressure  = 

1  Torr 

k  = 

1.6  X  lO'*  exp  (-18500/7) 

mole/m^s  Pa 

Twafer  = 

1,000  K 

Total  Flow  = 

0.2  slm 

kn  = 

0.19  Pa-®* 

Silane  fraction  = 

0.1 

h  = 

0.7  Pa* 

So,  this  simulation  would  require  conservation  equations  for  two  species  (silane  and  hydrogen) 
for  hydrogen  as  carrier  gas;  three  species  (silane,  H2  and  N2)  for  N2  as  carrier  gas. 

Figure  3  shows  flow  streamlines  in  the  reactor  for  H2  as  carrier  gas.  The  left  half  shows  a 
hypothetical  case  when  no  deposition  takes  place;  the  right  half  shows  the  actual  case  when 
deposition  is  in  progress.  In  both  cases,  the  general  flow  behavior  is  as  expected  from  inlet  to  the 
exit.  The  flow  is  laminar.  What  is  interesting  is  the  difference  between  the  two  cases  at  the 
substrate.  When  deposition  is  in  progress,  the  streamlines  are  normal  to  the  substrate,  indicating  a 
net  velocity  into  the  substrate,  which  is  due  to  a  mass  loss  following  surface  reactions  and 
consequent  deposition.  A  critical  inference  here,  as  pointed  out  in  [28],  is  that  the  temptation  to 
solve  for  the  fluid  flow  and  heat  transfer  first,  and  then  solve  for  species  conservation  using  the 
“frozen”  flow  field,  would  be  inappropriate;  it  would  give  wrong  velocity  profiles  near  the 
substrate  and  total  mass  flow.  The  species  conservation  is  essentially  coupled  to  fluid  flow. 

Figure  4  shows  a  similar  comparison  for  N2  as  carrier  gas.  For  the  same  mass  flow,  the  N2  has  a 
Reynolds  number  3.7  times  higher  than  that  for  the  H2  case.  Small  recirculation  zones  appear  in 
the  reactor.  Because  our  plotting  package  computes  the  streamlines  as  “particle  traces”,  the 
integration  error  makes  the  recirculation  zone  look  like  one  continuous  contour  line;  indeed,  it  is  a 
recirculation  region. 
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Figure  5  shows  the  temperature  contours  in  the  reactor,  comparing  the  effects  of  two 
types  of  carrier  gases.  As  expected,  there  is  a  strong  gradient  in  front  of  the  substrate.  The 
Reynolds  number  and  Peclet  number  for  the  N2  case  are  higher  than  for  the  H2  case.  As  a  result, 
conduction  spreads  out  the  temperature  a  little  farther  in  the  case  of  H2. 

Next,  we  discuss  species  mass  fraction  distribution  in  the  reactor.  Figure  6  shows  silane 
mass  fraction  in  a  silane/H2  mixture.  The  left  half  corresponds  to  a  case  which  does  not  include 
thermal  diffusion.  The  case  on  the  right  half  includes  thermal  diffusion.  Note  that  thermal 
diffusion  results  in  the  migration  of  species  due  to  temperature  gradients.  In  general,  if  molecular 
sizes  are  comparable,  the  heavier  species  concentrate  on  the  colder  regions  and  the  lighter  species 
migrate  to  warmer  regions.  In  these  reactors,  there  is  a  sharp  temperature  gradient  in  front  of  the 
susceptor,  as  seen  earlier.  As  a  result,  the  picture  of  silane  distribution  is  significantly  affected  by 
thermal  diffusion,  as  shown  in  Fig.  6.  Since  our  objective  is  to  validate  the  Ramshaw  model  for 
thermal  diffusion,  a  comparison  is  drawn  in  Fig.  7  for  Ramshaw  model  vs  empirical  coefficients. 
The  contours  are  qualitatively  similar.  The  quantitative  differences  are  actually  magnified  by 
surface  deposition  and  consumption/generation  of  species  at  the  substrate.  Indeed,  when  we 
artificially  turned  off  the  surface  deposition  boundary  condition,  the  comparison  of  the  two 
models,  in  terms  of  the  species  distribution,  was  good.  Figure  8  shows  the  growth  rate  of  silicon 
as  a  function  of  radius  for  all  of  the  above  three  cases.  The  Ramshaw  model  is  within  15%  of  the 
empirical  model.  The  agreement  is  much  closer  for  N2  as  carrier  gas.  Analogous  plots  are  shown 
in  Figs.  9  through  13  for  the  N2  case.  Note  that  this  is  a  three  component  system;  we  show  mass 
fractions  of  silane  and  H2  with  nitrogen  being  the  carrier  gas.  The  qualitative  behavior  and  the 
explanation  of  figures  are  the  same  as  above. 

The  above  investigation  reveals  that  the  model  proposed  by  Ramshaw  [12-14]  for  thermal 
diffusion  coefficients  is  reliable  and  can  be  used  in  CVD  studies.  We  also  used  this  opportunity  to 
validate  our  analysis  and  code  against  the  published  results  of  Kleijn,  et  al.  [28]  for  silicon.  After 
this  validation,  we  proceeded  to  analyze  CVD  of  high  temperature  materials  and  the  results  of  the 
investigations  are  presented  in  the  next  two  sections. 

6,2  CVD  of  silicon  carbide 

For  silicon  carbide  CVD,  many  different  source  gas  mixtures  can  be  used:  SiHVCHi, 
SfflU/CsHg,  CHsSi  CI3/H2,  methyldichlorosilane,  or  dimethyldichlorosilane.  In  general,  source 
gases  capable  of  delivering  C  and  Si  are  needed.  They  can  be  two  separate  feed  streams,  each 
containing  Si  or  C  individually,  as  in  silane-hydrocarbon  mixtures  or  chlorosilanes  which  contain 
both  Si  and  C.  A  detailed  discussion  on  SiC  preparation,  chemistry,  a  consolidated  historical 
review,  properties,  etc.,  can  be  found  in  Gmelin  Handbook  [29].  In  Phase  I,  we  compiled  reaction 
chemistry  for  all  of  the  above  feed  gases.  But,  due  to  the  limited  scope  of  Phase  I,  computations 
were  performed  for  only  one  system,  namely,  methyltrichlorosilane  (MTS)  /H2;  other  systems  will 
be  investigated  and  compared  in  Phase  II.  The  reaction  set  for  MTS/H2  system  is  given  in  Table 
II  [30-33]. 
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Table  IT;  Methyltrichlorosilane  (MTS)  reactions  for  SiC  deposition 
k  =  AT**  exp(-E/RT);  units:  cal,  mole,  cm,  s. 

1.0  elO  denotes  1x10*°.  M  is  third  body. 


Reaction 

A 

P 

E 

1. 

CbSiCHs 

<r> 

SiCl3  +  CH3 

7.6el4 

0 

69312.0 

2. 

SiCbH  +  H 

SiCl3  +  H2 

2.47el2 

0 

2543.5 

3. 

SiCbH 

SiCl2  +  HCl 

2.6ell 

0 

47000.0 

4. 

SiCL,  +  H 

SiCl3  +  HCl 

1.5el2 

0 

3400.0 

5. 

H  +  H  +  M 

H2  +  M 

2.95el8 

-1.0 

0 

6. 

CH4  +  H 

CH3  +  H2 

1.259el4 

0 

11900.0 

7. 

CH4  +  M 

<r> 

CH3  +  H  +  M 

1.413el7 

0 

88400.0 

8. 

CH3  +  CH3 

C2H6 

8.913el2 

0 

0 

9. 

C2H6  +  H 

C2H5  +  H2 

5.4e2 

3.5 

5210.0 

10. 

C2H5  +  M 

C2H4+H  +  M 

1.99el5 

0 

30000.0 

11. 

CH3+  CH3 

<r> 

C2H5  +  H 

2.8el3 

0 

13592.0 

12. 

C2H6  +  CH3 

C2H5  +  CH4 

5.5el 

4.0 

8300.0 

13. 

C2H4  +  H 

C2H5 

2.21el3 

0 

2066.0 

14. 

C2H4  +  M 

C2H2  +  H2  +  M 

1.5el5 

0 

55800.0 

Regarding  surface  processes,  at  present  we  account  for  them  using  a  sticking  coefficient 
approach,  though  we  have  initiated  a  more  detailed  procedure.  The  latter  involves  a  description 
of  all  elementary  steps  in  a  surface  activity  (such  as  species  adsorption,  surface  reaction,  product 
desorption,  etc.)  and  also  keeping  track  of  various  surface  sites  and  mass  balances  for  various 
adsorbates.  This  procedure  will  be  further  developed  and  fully  implemented  in  Phase  II. 

Typically,  all  radicals  are  very  active  on  the  surface  and  may  have  a  unity  sticking  coefficient.  The 
reactive  sticking  coefficients  for  the  more  stable  species  are  usually  smaller:  2  x  10'^,  2  x  10'^  and 
5  X  10'^  for  C2H4,  C2H2  and  CH4,  respectively. 

First,  we  describe  results  for  SiC  deposition  in  a  tubular  reactor.  The  motivation  comes 
from  commercial  processes  used  to  deposit  SiC  using  large  cylindrical  reactors  for  a  thickness  up 
to  1”  for  very  high  temperature  material  applications  [34,35].  On  a  much  smaller  scale,  CVI  can 
be  thought  of  as  CVD  in  a  cylindrical  pore.  In  any  case,  this  provides  us  an  opportunity  to 
compare  against  the  data  from  ref  [32  and  33].  The  cylindrical  reactor  is  30  cm  long  and  1.5  cm 
in  diameter.  The  operating  conditions  are  summarized  in  Table  III.  The  experiments  in  ref  [33] 
indicate  an  isothermal  middle  section  in  the  reactor. 
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Table  in.  SiC  deposition  in  a  tubular  reactor 


Pressure  =  1  atm 


Feed  gas;  methyltrichlorosilane/H2  =0.1 


Deposition  =  1,300  K 


Total  flow:  600  seem 


Figure  14  shows  the  concentration  of  hydrocarbon  species  along  the  length  of  the  reactor. 
Methane  is  the  most  abundant  hydrocarbon;  C2H4,  and  C2H2  are  the  other  stable  species.  C2H5 
concentration  is  too  small  to  show  on  the  scale  of  Fig.  14.  The  chlorosilane  concentrations  are 
shown  in  Fig.  15.  The  MTS  mole  fraction  decreases  fi'om  about  9%  at  the  inlet  to  5.5%  at  the 
exit.  Some  representative  reaction  rates  are  shown  in  Fig.  16.  Reaction  2  essentially  proceeds  in 
the  reverse  direction,  generating  SiCbH,  which  in  turn  decomposes  to  produce  SiCl2  radical.  The 
two  primary  radicals  from  MTS  decomposition  form  only  a  small  part  of  the  mixture,  due  to  their 
consumption  in  subsequent  reactions.  Temperatures  higher  than  1,300  K  would  favor  higher 
concentration  of  acetylene  than  obtained  here.  The  results  in  Figs.  14-16  are  in  good  agreement 
with  those  from  ref  [32,33]. 

Next  we  discuss  SiC  deposition  in  a  reactor  similar  to  the  one  considered  earlier  for 
silicon.  The  dimensions  of  the  reactor  are  shown  in  Fig.  17.  The  reactor  operating  conditions  are 
shown  in  Table  IV.  We  use  a  feed  gas  mixture  of  methyltrichlorosilane  and  hydrogen.  A 
parametric  study  varying  the  operating  pressure  (0.02,  0.05  and  0.1  atm)  was  conducted  to 
identify  its  influence  on  growth  rate  and  stoichiometry. 


Table  IV.  SiC  depostion 

Feed  gas:  MTS/H2  =  0.1 
Total  flow  rate  =  600  seem 


Temperature  =  1,300  K 

Pressure  =  0.02,  0.05  and  0.1  atm 


Flow  streamlines  in  the  reactor  at  0. 1  atm  are  shown  in  Fig.  18.  On  the  left,  we  show  a 
hypothetical  case  without  deposition  and  on  the  right,  the  actual  case  with  deposition.  In  both 
cases,  there  are  recirculation  zones.  As  pointed  out  earlier,  the  quirks  of  our  plotting  package 
make  the  recirculation  zone  look  like  one  continuous  line;  these  are  indeed  recirculation  zones. 
The  ones  near  the  wafer  are  due  to  buoyancy,  since  the  wafer  temperature  is  higher  than  the 
surroundings.  When  the  deposition  is  in  progress,  there  is  a  net  normal  velocity  at  the  wafer 
indicating  mass  loss.  As  discussed  earlier,  this  situation  precludes  the  “frozen  fluid  flow” 
approach  to  solve  the  species  equations.  Figure  19  shows  streamlines  at  0.02  and  0.05  atm,  with 
deposition  in  progress.  These  are  not  significantly  different  from  the  right  side  of  Fig.  18,  since 
the  mass  flow  in  all  three  cases  is  the  same.  The  temperature  contours  inside  the  reactor  at  two 
pressures  are  shown  in  Fig.  20.  Once  again,  qualitatively  they  are  very  similar. 
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Next,  we  inspect  mass  fraction  contours  of  some  key  reactive  species.  Mass  fraction  of 
MTS  at  two  pressures  is  shown  in  Fig.  21.  Note  that  the  decomposition  reaction  (reaction  1  in 
Table  II)  occurs  more  readily  near  the  susceptor,  due  to  local  high  temperatures.  The  minimum 
value  of  MTS  concentration  at  0. 1  atm  is  smaller  than  that  at  lower  pressures.  This  is  due  to  the 
fact  that  the  residence  time  is  smaller  as  the  pressure  is  decreased,  following  an  increase  in 
velocity.  As  such,  MTS  does  not  decompose  completely  as  the  pressure  is  reduced.  Mass 
fraction  contours  for  SiCl2  and  CHU  are  shown  in  Figs.  22  and  23,  respectively.  Again,  the 
fractions  of  these  species  are  higher  at  higher  pressures.  Their  concentrations  decrease  rapidly  as 
one  moves  away  from  the  susceptor.  Consistent  with  the  above  species  distributions,  the  growth 
rate  increases,  with  an  increase  in  pressure,  as  shown  in  Fig.  24.  The  rates  and  dependence  on 
pressure  shown  here  are  similar  to  those  reported  in  [32,33].  The  peak  in  the  rates  in  Fig.  24 
coincides  with  the  radial  location,  where  peaks  in  species  concentrations  are  seen.  Finally,  the 
deposit  stoichiometry  is  shown  in  Fig.  25  as  a  function  of  pressure.  The  Si  fraction  decreases 
with  pressure,  as  seen  in  experiments  [33];  however,  we  underpredict  here  the  Si  atomic  fraction 
due  to  the  values  of  sticking  coefficients  used.  A  more  detailed  surface  chemistry  planned  for 
Phase  II  would  improve  this  prediction. 

6.3  Boron  deposition  on  a  tungsten  wire 

Boron  is  a  hard,  light,  refractory  material  that  exhibits  high  corrosion  resistance. 
Consequently,  there  are  many  applications  in  refractory  coatings  and  composite  materials.  One  of 
the  common  processes  is  boron  deposition  on  wires  of  tungsten,  tantalum,  molybdenum,  or 
titanium.  Here  we  consider  deposition  of  boron  on  tungsten  wires.  The  reactor  consists  of  a  30 
cm  long  Pyrex  tube  whose  radius  is  3  cm.  The  tungsten  wire  radius  is  0.01  cm.  The  reactor  can 
be  operated  in  a  batch  or  continuous  mode,  i.e.,  the  wire  may  be  stationary  or  moving.  Here,  the 
reactor  is  vertical.  The  operating  conditions  are  given  in  Table  V. 


Table  V;  Boron  deposition  parameters 
Pressure  =  50  K  Pa  BCI3/H2  mixture 

Wire  Temperature  =  1,000  K  BCI3  =  5% 

Total  flow  =  1,700  seem 


Normal  operation  of  the  reactor  requires  ohmic  or  resistive  heating  of  the  tungsten  wire  by 
applying  an  electric  current.  The  wire  transmits  heat  to  the  flowing  gas  by  conduction  and 
convection  and  also  loses  heat  to  the  surroundings  by  radiation.  The  ability  to  absorb  power  by 
the  wire  is  determined  by  the  resistivity  of  the  wire,  which  can  constantly  change  as  deposit  builds 
up.  All  of  this  physical  phenomena  have  been  accounted  for  in  the  wire  boundary  condition.  The 
boron  deposition  is  treated  using  a  surface  reaction; 
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BCb  +  3/2  H2  ->B  +  3HC1. 


The  rate  constant  is  given  by  2.0  exp  (-7594/T)  cm/s  [36-38], 

The  flow  inside  the  reactor  is  laminar.  The  thermal  contours  are  shown  in  Fig.  26.  Note 
that  the  computation  is  two-dimensional  axisymmetric.  The  gas  enters  at  the  bottom  at  room 
temperature  and  is  heated  progressively  as  it  reaches  the  top.  The  reactor  wall  temperature  also 
increases,  since  it  is  determined  by  the  radiation  it  receives  from  the  wire  and  heat  loss  to  the 
surroundings  by  convection  and  radiation.  Figure  27  shows  mass  fraction  contours  of  BCI3  and 
HCl;  Fig.  28  displays  H2  contours.  Only  a  small  amount  of  BCI3  is  lost  due  to  decomposition. 
There  is  also  a  gradient  in  BCI3  in  the  radial  direction.  HCl  formation  near  the  inlet  is  negligible, 
but  its  composition  increases  as  the  gas  moves  up  the  reactor.  As  expected,  its  concentration  is 
higher  near  the  wire  than  in  the  reactor  wail  region.  The  growth  rate  of  boron  was  found  to  be 
nearly  uniform  (except  for  the  inlet  region)  and  at  2.84  nm/s.  These  computations  demonstrate 
our  capability  to  model  moving  wire  deposition  problems. 
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7.  CVI  Results 


7.1  Reactor  Geometry 

In  Phase  I,  CVI  of  a  slab-like  preform  using  microwave  heating  is  studied;  we  used  the 
Phase  I  opportunity  to  add  to  our  code  a  module  to  solve  microwave  heating  equations. 
Conventional  heating  does  not  require  this  effort,  but  needs  only  an  energy  conservation  equation. 
The  new  capability  is  demonstrated  here  using  the  following  case  study.  Phase  II  will  involve 
studies  of  infiltration  of  large  objects,  two-dimensional  analysis  of  the  preform,  and  conventional 
fijmace  heating  in  isothermal  CVI. 

The  CVI  system  investigated  is  shown  in  Figure  29.  It  consists  of  a  free-standing  preform 
of  slab  geometry  composed  of  cylindrical  fibers  oriented  randomly  in  3-dimensional  space.  The 
preform  resides  in  a  cavity  where  microwaves  are  incident  onto  both  faces  of  the  preform.  A 
gaseous  mixture  of  known  composition  and  temperature  flows  continuously  through  the  cavity, 
keeping  a  constant  gas  composition  at  the  preform  surface.  Reactants  infiltrate  the  fibrous 
structure  and  participate  in  a  reaction  network,  resulting  in  the  deposition  of  a  solid  material  that 
fills  the  pores.  Deposition  continues  until  complete  densification  is  obtained  or  stops  if  premature 
pore  closure  occurs  at  the  peripheries.  The  goals  are  as  follows:  for  arbitrary  (but  reasonable) 
chemistry,  determine  the  spatiotemporal  profiles  of  composite  mass  density,  temperature,  species 
concentration  and  pressure  distribution.  An  important  optimization  question  is  the  following: 
what  is  the  optimum  heating  microwave  power,  or  optimum  heating  schedule  (by  varying  power 
as  a  function  of  time)  that  results  in  complete  densification  in  the  least  amount  of  time? 

7.2  Chemical  Kinetics 


The  deposition  of  silicon  carbide  by  decomposition  of  methyltrichlorosilane  (MTS)  was 
taken  as  the  chemical  system  for  investigation.  Although  this  system  has  been  the  subject  of 
several  works,  most  studies  have  assumed  (due  to  the  lack  of  kinetic  data)  an  overall  deposition 
reaction  that  is  first  order  in  MTS  concentration  [39].  Recently,  however,  several  workers  [32, 

40,  41]  have  investigated  in  some  detail  the  gas  phase,  as  well  as  the  surface  reactions  leading  to 
the  deposition  of  SiC.  These  studies  have  focused  on  the  homogeneous  [32],  as  well  as 
heterogeneous  processes  taking  place  [41].  Based  on  the  experimental  investigation  conducted  by 
Tsai,  et  al.  [41],  the  following  mechanism  was  used  to  describe  the  deposition  process. 


CHjSiCls  +  H2  ^SiCl2  +  CH4  +  HCl 
SiCb  +  *  ^  SiClz* 

CH4  +  *  ->  CH4* 

SiClj*  +CH4*  ^  SiC  +  H2  +  2HC1* 

HCl*  ->  HCl  +  * 

In  the  above  reactions,  *  denotes  a  site  on  the  surface.  A  species  with  an  *  denotes  an  adsorbed 
species.  This  mechanism  postulates  that  the  precursor  decomposes  in  the  gas  phase  to  form 
intermediate  products  which  then  adsorb  on  the  surface.  The  adsorbates  react  on  the  surface  to 
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yield  SiC  deposit  and  hydrogen  and  hydrogen  chloride  by-products.  We  realize  that  this  is  a 
greatly  simplified  picture  of  the  actual  chemical  mechanisms  taking  place.  However,  we  feel  that 
this  mechanism  captures  the  salient  features  (decomposition,  adsorption,  surface  reaction, 
desorption)  of  a  variety  of  different  chemical  precursor  systems.  In  the  absence  of  reliable  kinetic 
data,  any  more  complicated  mechanism  is  not  warranted  at  the  present  time.  However,  we  would 
like  to  stress  that  the  model  formulation  is  now  general  enough  to  accommodate  any  number  of 
chemical  species  and  any  type  of  reactions.  Hence,  the  model  can  be  used  for  any  CVI  systems 
with  known  chemistry. 

For  the  reaction  mechanism  given  above,  a  standard  Langmuir-Hinshelwood  analysis 
yields  the  following  expressions  for  the  rate  of  MTS  decomposition  and  the  rate  of  SiC 
deposition,  respectively. 


Rl  —  ki  Ch2CmtS  -  k-i  CsiCI2  CcH4  CHI 
R2  =  k2  CcH4  CsiC12/(l  +  kal  CsiC12  +  ka2  CcH4)^ 

The  corresponding  rate  coefficients  have  been  measured  by  Tsai,  et  al.,  [41]: 

k,  =  2.0  X  lO'^^expC-  448.2/RT)  m^/(mol  s) 
k-i  =  1.1  X  10^°exp(-416.2/RT)  mV(mol  s) 
k2-  1.633  X  10'*exp(-319.8/RT)  m\mol  s) 
kai  =  5.0  X  10‘*exp(-21.6/RT)  mV(mol) 
ka2  =  7.1  X  10^exp(-33.1/RT)  mV(mol) 

In  accordance  with  the  reaction  network  above,  five  gaseous  components  were  included  in  the 
model;  namely,  CHsSiCb,  H2,  SiCU,  CH4,  and  HCl.  The  reaction  network  can  be  augmented  to 
include  etching  of  the  deposit  by  HCl,  which  can  become  important  at  high  temperatures. 

7.3  Results 


In  the  following,  attention  is  focused  on  the  influence  of  microwave  power  on  the 
spatiotemporal  behavior  of  the  densification  process.  Base  parameter  values  used  for  calculations 
are  shown  in  Table  VI. 
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Table  VI  Base  parameter  values  used  for  calculation 


Symbol 

Name 

Base  Value 

L 

Half  preform  thickness 

1  mm 

‘|)a,o 

Initial  accessible  porosity 

0.5 

h 

heat  transfer  coefficient 

329  WWK) 

Po 

Initial  pressure 

1  atm 

ff 

Fiber  radius 

4.0  pm 

Tb 

Ambient  temperature 

300  K 

r 

Reference  temperature 

300  K 

Xib 

CHsSiCls  mole  fi’action 

0.1 

X2b 

Hydrogen  mole  fraction 

0.9 

lo 

Incident  power  flux 

1.5MW/m^ 

f 

Frequency 

2.45  GH^ 

k’ 

Relative  dielectric  constant 

60 

k” 

Relative  loss  factor 

35 

k^ 

reference  thermal  conductivity 

6.58  X  10’^ 

J/(m  s  K) 

E 

Young’s  modulus  of  elasticity 

400  GPa 

V 

Poisson’s  ratio 

0.2 

ttT 

Coefficient  thermal  expansion 

5.0x  lO-^K'* 

Figure  30  shows  the  transient  behavior  of  the  porosity  profiles  for  an  incident  power  flux 
of  1.5  MW/m^.  One  can  see  that  the  porosity  is  quite  uniform  at  the  initial  stages  of  the  process. 
As  deposition  continues,  however,  the  temperature  increases  significantly  (Figure  31)  due  to 
further  microwave  absorption  by  the  denser  composite.  This  temperature  increase  causes  faster 
reaction  rate  and  reactant  depletion  as  the  reactants  diffiise  into  the  preform.  Eventually,  sealing 
of  the  preform  surface  results  in  residual  accessible  porosity  trapped  inside  near  the  composite 
center.  The  system  behavior  reflects  the  strong  interplay  between  microwave  absorption  and  SiC 
deposition  and  its  effect  on  density  uniformity.  Note  that  the  conditions  of  Fig.  30  lead  to  very 
rapid  deposition;  about  half  of  the  preform  is  filled  within  10  mins! 

Reducing  the  incident  power  flux  (Fig.  32)  can  alleviate  reactant  depletion  problems 
because  the  composite  temperature  is  inevitably  lower.  Under  these  conditions,  the  deposition 
process  follows  an  “inside-out”  desensification  pattern  up  to  about  1  hr.  Thereafter,  the 
microwave  absorption  feedback  mechanism  results  once  again  in  high  enough  temperature  to 
cause  premature  surface  pore  closure.  Investigation  of  still  lower  incident  power  fluxes  revealed 
an  ignition-type  behavior  characterized  by  almost  no  deposition  taking  place  below  a  certain 
incident  flux  (about  1.15  MW/m^).  Above  this  value,  however,  the  system  temperature  increases 
slowly  at  early  times,  but  then  reaches  a  point  of  accelerated  rate  of  increase. 
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One  may  consider  lower  processing  pressures,  as  shown  in  Figure  33  (Figs.  30-32  were 
for  a  pressure  of  1  atm).  Although  the  deposit  uniformity  has  improved  at  the  expense  of  longer 
processing  time,  entrapment  of  porosity  still  occurs  due  to  more  microwave  power  absorption  as 
densification  progresses 

All  the  above  examples  were  for  a  constant  incident  microwave  power.  The  picture  can 
be  improved  by  using  a  variable  power  schedule.  For  example,  one  can  use  a  high  power  initially 
to  fill  in  most  of  the  preform  and  then,  at  some  point  during  the  infiltration,  switch  to  a  lower 
power  to  avoid  premature  pore  closure  at  the  surface  of  the  preform.  An  example  of  such  a  step 
change  in  power  is  shown  in  Figure  34.  The  deposit  uniformity  is  substantially  improved  and  an 
inside-out  densification  pattern  is  obtained.  However,  it  may  be  difficult  to  determine  a  priori  the 
proper  switching  time(s). 

An  alternative  power  scheduling  technique  that  has  been  found  to  lead  to  an  increase  in 
density  uniformity  without  compromising  the  processing  time  is  pulsed  microwave  power  [21].  In 
this  technique,  the  source  power  is  modulated  in  time  with  a  specific  period  and  duty  cycle. 

During  the  low-power  part  of  the  cycle,  the  temperature  of  the  composite  decreases,  reducing  the 
reaction  rate  and  thus  allowing  the  reactants  to  penetrate  into  the  composite.  This  alleviates 
diffiisional  limitations  within  the  composite,  minimizing  density  nonuniformities.  The  high-power 
part  of  the  cycle  leads  to  rapid  reaction  rates,  thereby  minimizing  the  processing  time. 

Figure  35  illustrates  the  time  history  of  the  spatially-averaged  accessible  porosity  (right  y- 
axis),  as  well  as  that  of  the  center  (Tc),  and  surface  (Ts)  temperature  of  the  composite  for  a  case 
with  square  wave  modulation  of  the  power.  The  pulse  period  was  180  s  with  a  duty  cycle  of 
50%.  The  high  and  low  power  levels  were  1.3  and  1.2  MW/m^,  respectively.  During  the  early 
stages  of  the  process  the  preform  temperature  increases  from  its  initial  value  of  300  K  to  about 
1200  K  without  any  substantial  deposition  of  SiC  taking  place.  As  densification  progresses,  the 
composite  temperature  keeps  increasing  with  time,  but  does  not  reach  the  point  of  inducing 
surface  power  closure.  Power  modulation  yields  an  “inside-out”  densification  pattern,  as  shown 
in  Fig.  36.  Note  that  the  processing  time  is  rather  short  (hours)  compared  to  that  of  Fig.  34  (that 
also  resulted  in  inside-out  deposition). 

The  gaseous  species  mole  fraction  profiles  in  the  composite  for  the  above  power- 
modulation  conditions  are  depicted  in  Figure  37  at  two  time  instants.  It  can  be  seen  that  MTS 
decomposes  to  a  large  extent  during  the  high-power  part  of  the  cycle  (at  t  =  6900  s,  dashed  line), 
resulting  in  the  formation  of  SiCl2  and  CHt.  These  species  subsequently  react  heterogeneously  to 
form  SiC.  During  the  low-power  part  of  the  cycle  (t  =  6600  s,  solid  line)  the  temperature 
decreases  and  hence,  MTS  is  allowed  to  infiltrate  deeper  into  the  composite.  This  results  in  better 
deposit  uniformity.  It  should  be  noted  that  substantial  amounts  of  HCl  are  accumulated  within  the 
composite.  It  is  conceivable  that  HCl  will  etch  back  some  of  the  deposit.  It  is  straightforward  to 
incorporate  an  etching  reaction  in  the  model,  which  will  be  done  in  Phase  II. 

The  above  computations  demonstrate  how  one  can  use  numerical  simulations  to  conduct 
design  studies.  Here,  due  to  the  limited  scope  of  Phase  I,  we  restricted  our  attention  only  to  the 
heating  mode  as  a  design  variable.  In  Phase  II,  all  other  possible  variables,  such  as  size  and  nature 
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of  preform,  infiltration  pressure  and  temperature  and  their  gradients,  if  any,  precursor  feed  rates, 
external  reactor  conditions,  will  be  considered.  Also,  other  systems  of  interest,  such  as  C-C,  will 
be  included.  A  systematic  comparison  of  different  available  source  gases  for  a  system  can  also  be 
made. 
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8.  Concluding  Remarks  and  Directions  for  Phase  n 


The  objective  of  this  multiphase,  multiyear  STTR  project  is  to  develop  computational 
simulation  and  codes  for  structural  materials  and  composites  processing  which  would  aid  in 
understanding  of  mechanisms,  reactor  development,  process  optimization,  and  process  control, 
and  serve  as  a  complementary  design  tool.  There  are  a  number  of  techniques  used  in  structural 
materials  processing:  hot  pressing,  sintering,  thermal  or  plasma  spraying,  CVD,  etc.  We  chose 
here  two  of  the  most  widely  used  techniques:  CVD  and  CVI,  which  are  closely  related  in  many 
ways.  Also,  the  modeling  focus  here  is  on  what  is  happening  inside  the  reactor — physical  and 
chemical — with  a  view  to  correlate  process  variables  (for  which  there  are  knobs  on  the  reactor 
panel,  such  as  pressure,  flow  rates,  heating/cooling,  etc.)  to  process  performance  figures-of-merit, 
such  as  growth  rate,  uniformity,  surface  composition,  etc.  This  exercise,  in  itself,  is  significantly  a 
large  task,  with  much  value  to  the  process  engineer.  While  it  is  desirable  to  extend  the  analysis  to 
the  microtopographical  evolution  of  the  coating  or  film,  etc.,  in  CVD  and  also  model  the  bulk 
mechanical  and  thermal  properties,  we  must  realize  that  it  would  be  impossible  to  finish  all  of  the 
tasks  in  one  set  of  Phase  I  and  II  programs.  The  latter  is  a  large  field  in  itself,  and  currently 
addressed  by  a  number  of  researchers  in  the  community.  On  the  other  hand,  we  did  consider  pore 
or  microtopographical  evolution  in  CVI,  which  otherwise  would  be  incomplete. 

The  Phase  I  work  is  to  demonstrate  the  concept  that  modeling  and  simulation  can  provide 
insight  into  processing  and  help  process  design.  This  has  been  successfully  demonstrated.  The 
accomplishments  are  as  follows: 

*  An  existing  SRA  computational  fluid  dynamics  code,  well-tested  for  a  variety  of  flow 
problems,  has  been  adapted  for  the  study  of  CVD  and  CVI  processes.  Multicomponent 
species  transport  for  an  arbitrarily  high  number  of  species,  thermal  diffusion  effects, 
homogeneous  and  heterogeneous  reactions,  proper  deposition  boundary  conditions, 
analysis  of  heat  transfer  within  solid  structures,  and  other  features  peculiar  to  chemically 
reacting  flows  have  been  included. 

*  Lower  order  models  to  analyze  a  complex  set  of  chemical  reactions — involving  hundreds 
of  species  in  hundreds  of  reactions — have  been  further  developed.  They  were  used  to 
conduct  a  sensitivity  study  of  available  reaction  data  before  using  in  the  multidimensional 
simulation. 

*  A  module  to  analyze  chemical  vapor  infiltration  (CVI)  process  has  been  developed  by  the 
STTR  university  partner. 

*  The  codes  and  capabilities  were  demonstrated  in  several  preliminary  case  studies:  CVD  of 
silicon  carbide,  CVD  of  boron  on  fibers,  and  CVI  of  SiC  in  slab  preforms.  In  each  case, 
the  model  has  provided  insight  into  transport  and  reaction  mechanisms  of  the  processes 
which  is  essential  for  process  optimization.  For  example,  in  SiC  deposition,  the  effect  of 
operating  pressure  on  deposition  rate  and  SiC  stoichiometry  was  studied.  In  the  case  of 
SiC-SiC  CVI,  the  evolution  of  the  preform  infiltration  characteristics  was  studied,  as  a 
function  of  the  heating  profile,  in  an  attempt  to  optimize  heating  schedules  to  avoid 
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premature  pore  closing.  These  are  some  preliminary  examples  of  how  a  well-tuned  model 
can  be  used  as  an  adjunct  to  experimental  efforts  to  increase  our  understanding  of  process 
mechanisms.  Comparison  against  experimental  results,  wherever  possible,  showed  that 
the  model  reproduces  observed  behavior. 

The  success  of  the  Phase  I  effort  is  encouraging.  In  line  with  the  original  STTR 
Solicitation,  SRA  would  further  develop  the  model  and  codes  to  simulate  CVD  and  CVI  in 
structural  materials  processing.  Phase  II  tasks  would  include:  completion  of  CVD  model 
development,  improvement  of  numerical  procedures  for  robustness  and  rapid  convergence, 
completion  of  CVT  module  and  two-dimensional  description  of  transport  within  preform,  coupling 
of  reactor  and  CVI  modules  for  complete  CVI  analysis,  demonstration  of  several  processes  of 
interest  to  the  Air  Force  Laboratories  and  aerospace  industry,  and  a  thorough  model  validation  by 
comparison  against  experiments.  We  have  made  arrangements  with  three  large  companies  to 
provide  us  with  experimental  data.. 

The  Phase  III  effort  would  aim  at  developing  this  code  into  a  marketable  product.  After 
consulting  the  companies  and  DoD  personnel  we  have  interacted  with  in  the  past  year,  the 
following  strategy  has  been  developed  for  marketing.  There  are  many  fluid  dynamics  softwares 
available  in  the  market  today.  They  are,  however,  general  purpose  codes;  for  example,  the  same 
code  is  marketed  to  chemical,  semiconductor,  and  aerospace  industries  to  solve  generic  Navier- 
Stokes  equation-based  problems.  Though  these  codes  come  with  a  user-interface,  they  can  now 
be  used  only  by  experts  in  simulation.  Many  small  companies  and  even  large  corporations  no 
longer  have  staff  dedicated  to  modeling.  If  modeling  were  to  help  industry,  the  task  nowadays 
needs  to  be  done  by  process  engineers.  Often  process  engineers  are  not  experts  in  modeling  and 
simulation;  they  lack  the  background  in  several  physical  phenomena,  such  as  theories  behind 
multicomponent  diffusion,  thermal  diffusion,  reaction  mechanisms,  surface  processes,  etc.  Worse, 
they  do  not  have  the  time  and  resources  to  search  for  the  enormous  amount  of  transport  and 
kinetic  data  essential  for  the  process  in  hand.  So  we  learned  from  a  survey  of  industrial 
colleagues  that  we  must  develop  the  above  data  for  key  processes  of  interest  to  the  materials 
community  and  make  it  part  of  the  interface.  After  all,  there  are  only  a  handful  of  “proven  and 
new”  materials  and  chemistries  for  which  we  need  to  assemble  every  bit  of  information  that  is  an 
integral  part  of  the  simulation  input.  Also,  the  code,  solution  techniques,  and  the  user-interface 
would  be  specific  to  the  materials  problems.  This,  we  were  told,  would  make  it  attractive  to 
process  engineers  and  allow  SRA  to  sell  or  lease  the  software. 
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ier  gas.  Small  recirculation 
facts  make  the  recirculation  z 


Fig.  5  Temperature  contours.  Left:  N  carrier  gas.  Right:  H  carrier 
gas.  Reactor  conditions  as  in  Table  I.  Contour  interval  =  50  K 
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Growth  rate  as  a  function  of  radius  for  H2  as  carrier  gas 
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14  Hydrocarbon  species  concentration  in  a  tubular  reactor. 
SiC  deposition  conditions  as  in  Table  III. 


Fig.  15  Chlorosilane  species  concentration  in  a  tubular  reactor. 
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Fig,  16  Reaction  rates  of  some  key  reactions  in  Table  II.  Operating 
conditions  are  as  in  Table  III. 


Fig.  17  Schematic  of  the  SiC  deposition  reactor.  Dg  =  5  cm.  Dp  =  1  cm 
and'  Zp  =  0.5  cm.  Susceptor  temperature  =  1300  K. 


Fig,  18  Flow  streamlines  in  a  SiC  deposition  reactor  at  p  =  0,1  atm.  Left:  case 
without  deposition.  Right:  with  deposition.  Note  the  recirculation  in 
the  reactor;  the  plotting  program  makes  them  look  like  one  contour  line, 
but  it  is  recirculation. 


Fig,  19  Flow  streamlines  in  SiC  depostion.  Left:  0.05  atm.  Right:  0.02  atm. 
See  comments  on  Fig.  18. 
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Fig,  26  Temperature  contours  in  a  boron  deposition  reactor 
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Fig. 27  Mass  fraction  contours:  Boron  deposition 
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Fig.  28  Mass  fractidn  contours:  Boron  deposition 


Schematic  of  the  preform  in  SiC-SiC  CVI 


Fig.  30  Porosity  evolution.  Power  flux  =  1.5  MW/M^. 
Other  parameters  as  in  Table  VI. 
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Fig.  31  Temperature  contours  corresponding  to  Fig.  30. 


Fig.  32  Porosity  evolution  for  a  power  flux  of  1.3 
Parameters  as  in  Table  VI. 
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Fig.  33  Porosity  evolution  at  reduced  pressure  (0.5  atm). 
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Fig*  34  Effect  of  a  step  change  in  applied  power  on  porosity  evolution 
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Fig.  35  Variation  of  the  spatially-averaged  porosity,  centerline 


temperature  and  surface  temperature  when  the  input  power 
has  a  square-wave  modulation.  All  other  parameters  as  in 
Table  VI. 
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Fig.  36  Porosity  evolution  corresponding  to  Fig.  35 
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Fig.  37  Species  mole  fraction  at  two  time  instants  for  conditions 
in  Fig.  35. 


